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INTRODUCTION 


This  Appendix  contains  the  reprints  published  under  JSEP  in  the  time  September  1993 
to  September  1994. 

In  addition  to  the  8  reprints  contained  herein,  there  are  10  papers  already  accepted 
for  publication  during  the  next  contract  period,  12  papers  submitted  and  10  papers  in 
preparation. 
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JSEP  REFEREED  JOURNAL  PAPERS 
PUBLISHED  SEPTEMBER  1993  TO  SEPTEMBER  1994 

1.  G.A.  Somers  and  P.H.  Pathak,  “Efficient  Numerical  and  Closed  Form  Asymptotic 
Representations  of  the  Dyadic  Aperture  Green’s  Function  for  Material  Coated  Ground 
Planes,”  Radio  Science,  vol.  29,  pp.  465-481,  March-April  1994. 

2.  H.C.  Ly  and  R.G.  Rojas,  “EM  Plane  Wave  Diffraction  by  a  Material  Coated  Per¬ 
fectly  Conducting  Half-Plane — Oblique  Incidence,”  Microwave  and  Optical  Technology 
Letters,  vol.  7,  no.  6,  pp.  262-266,  April  1994. 

3.  R.  Torres  and  E.H.  Newman,  “Integral  Equation  Analysis  of  a  Sheet  Impedance  Coated 
Window  Slot  Antenna,”  IEEE  Transactions  on  Antennas  and  Propagation,  vol.  42, 
pp.  541-544,  April  1994. 

4.  J.L.  Blanchard,  E.H.  Newman  and  M.E.  Peters,  “Integral  Equation  Analysis  of  Arti¬ 
ficial  Media,”  IEEE  Transactions  on  Antennas  and  Propagation,  vol.  42,  pp.  727-731, 
May  1994. 

5.  R.  Lee  and  V.  Chupongstimun,  “A  Partitioning  Technique  for  the  Finite  Element 
Solution  of  Electromagnetic  Scattering  from  Electrically  Large  Dielectric  Cylinders,” 
IEEE  Transactions  on  Antennas  and  Propagation,  vol.  42,  pp.  737-741,  May  1994. 

6.  R.  Lee  and  T.T.  Chia,  “Analysis  of  Electromagnetic  Scattering  from  a  Cavity  with  a 
Complex  Termination  by  Means  of  a  Hybrid  Ray-FDTD  Method,”  IEEE  Transactions 
on  Antennas  and  Propagation,  vol.  41,  pp.  1560-1569,  November  1993. 

7.  J.O.  Jevtic  and  R.  Lee,  “On  the  Choice  of  Metrons  in  the  MEI  Method,”  COMPEL  — 
Int.  Journal  Comp,  and  Math,  in  Electrical  and  Electronic  Eng.,  vol.  13,  Supplement 
A,  pp.  217-222,  1994. 

8.  J.O.  Jevtic  and  R.  Lee,  “A  Theoretical  and  Numerical  Analysis  of  the  Measured 
Equation  of  Invariance,”  IEEE  Transactions  on  Antennas  and  Propagation,  vol.  42, 
pp.  1097-1105,  August  1994. 
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JSEP  RELATED  REFEREED  JOURNAL  PAPERS 
ACCEPTED  FOR  PUBLICATION 
SEPTEMBER  1993  TO  SEPTEMBER  1994 


1.  U.  Pekel  and  R.  Lee,  “An  A  Posteriori  Error  Reduction  Scheme  for  the  Three  Di¬ 
mensional  Finite  Element  Solution  of  Maxwell’s  Equations,”  IEEE  Trans.  Microwave 
Theory  and  Techniques. 

2.  M.E.  Peters  and  E.H.  Newman,  “Method  of  Moments  Analysis  of  Anisotropic  Artifi¬ 
cial  Media  Composed  of  Dielectric  Wire  Objects,”  IEEE  Transactions  on  Microwave 
Theory  and  Techniques  . 

3.  L.M.  Chou,  R.G.  Rojas  and  P.H.  Pathak,  “WH/GSMT  Based  Full- Wave  Analysis  of 
Planar  Transmission  Lines  Embedded  in  Mtiltilayered  Dielectric  Substrates,”  IEEE 
Trans,  on  Microwave  Theory  and  Techniques. 

4.  R.G.  Rojas,  “Integral  Equations  for  the  EM  Scattering  by  Homogeneous/Inhomogeneous 
Two  Dimensional  Chiral  Bodies,”  lEE  Proceedings-H. 

5.  M.F.  Otero  and  R.G.  Rojas,  “Synthesis  of  the  Frequency  Response  of  an  Inhomoge¬ 
neous  Resistive  Strip,”  Annales  of  Telecommunications^  Special  Issue  on  Radar  Cross 
Section. 

6.  H.T.  Anastassiu  and  P.H.  Pathak,  “High  Frequency  Analysis  of  Gaussian  Beam  Scat¬ 
tering  by  a  Two-Dimensional  Parabolic  Contour  of  Finite  Width,  ”  Radio  Science. 

7.  P.  Munk  and  P.H.  Pathak,  “EM  Scattering  by  a  Dielectric  Filled  Rectangular  Antenna 
Cavity  Recessed  in  a  Ground  Plane  and  Backed  with  an  Array  of  Loaded  Dominant 
Mode  Waveguides,”  IEEE  Transactions  on  Antennas  and  Propagation. 

8.  P.  Munk  and  P.H.  Pathak,  “A  Useful  Approximate  Analysis  of  the  EM  Scattering 
by  a  Rectangular  Antenna  Cavity  containing  an  Array  of  Dominant  Mode  Waveguide 
Loaded  Slots,”  Journal  of  EM  Waves  and  Applications,  (special  issue  on  EM  Scatter¬ 
ing). 

9.  P.  Rousseau  and  R.J.  Burkholder,  “A  Hybrid  Approach  for  Calculating  the  Scattering 
from  Obstacles  within  Large  Open  Cavities,”  IEEE  Transactions  on  Antennas  and 
Propagation. 

10.  W.P.  PineUo,  R.  Lee  and  A.C.  Cangellaris,  “Finite  Element  Modeling  of  Electromag¬ 
netic  Wave  Interactions  with  Periodic  Structures,”  IEEE  Trans.  Microwave  Theory 
and  Techniques. 
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JSEP  RELATED  PAPERS 
SUBMITTED  FOR  PUBLICATION 
SEPTEMBER  1993  TO  SEPTEMBER  1994 

1.  T.T.  Chia,  R.J.  Burkholder  and  R.  Lee,  “The  Application  of  FDTD  in  Hybrid  Methods 
for  Cavity  Scattering  Analysis,”  submitted  to  IEEE  Trans.  Antennas  and  Propagat. 

2.  T.L.  Barkdoll  and  R.  Lee,  “Finite  Element  Analysis  of  Bodies  of  Revolution  using  the 
Measured  Equation  of  Invariance,”  submitted  to  Radio  Science. 

3.  J.O.  Jevtic  and  R.  Lee,  “The  Relationship  Between  the  Reflection  Coefficient  and  the 
Error  for  Absorbing  Boundary  Conditions,”  submitted  to  Electromagnetics. 

4.  J.O.  Jevtic  and  R.  Lee,  “An  Investigation  of  the  Convergence  Properties  of  the  Mea¬ 
sured  Equation  of  Invariance,”  submitted  to  IEEE  Transactions  on  Antennas  and 
Propagation. 

5.  J.O.  Jevtic  and  R.  Lee,  “How  Invariant  is  the  Measured  Equation  of  Invariance?” 
submitted  to  IEEE  Microwave  and  Guided  Wave  Letters. 

6.  Y.S.  Choi-Grogan,  K.  Eswar,  P.  Sadayappan  and  R.  Lee,  “Sequential  and  Parallel  Im¬ 
plementations  of  a  Partitioning  Finite  Element  Method,”  submitted  to  IEEE  Trans¬ 
actions  on  Antennas  and  Propagation. 

7.  U.  Pekel  and  R.  Lee,  “A  3-D  Application  of  the  Bymoment  Method  for  Electromagnetic 
Scattering,”  submitted  to  Journal  of  Electromagnetic  Waves  and  Application. 

8.  E.H.  Newman  and  M.  Kragalott,  “Moment  Method  Analysis  of  the  Electric  Shielding 
Factor  of  a  Conducting  TM  Shield  at  ELF,”  submitted  to  IEEE  Trans,  on  Electro¬ 
magnetic  Compatability. 

9.  L.M.  Chou  and  R.G.  Rojas,  “Dispersion  and  Lateral  Leakage  of  Conductor  Backed 
Coplanar  Waveguides  with  Layered  Substrate  and  Finite-Extent  Lateral  Ground 
Planes,”  submitted  to  IEEE  Transactions  on  Microwave  Theory  and  Techniques. 

10.  P.  Rousseau  and  P.H.  Pathak,  “Time  Domain  Uniform  Geometrical  Theory  of  Diffrac¬ 
tion  (TD-UTD)  for  a  Curved  Wedge,”  submitted  to  IEEE  Transactions  on  Antennas 
and  Propagation. 

11.  F.  Obelleiro,  J.L.  Rodriguez  and  R.J.  Burkholder,  “An  Iterative  Physical  Optics  Ap¬ 
proach  for  Analyzing  the  Electromagnetic  Scattering  by  Large  Open-Ended  Cavities,” 
submitted  to  IEEE  Transactions  on  Antennas  and  Propagation. 

12.  P.H.  Pathak,  H.T.  Chou  and  R.J.  Burkholder,  “Ray  Like  Gaussian  Basis  Functions  for 
Analyzing  Propagation  into  and  Scattering  from  Large  Open  Cavities  with  the  GRE.” 
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JSEP  RELATED  PAPERS 
IN  PREPARATION  FOR  PUBLICATION 
SEPTEMBER  1993  TO  SEPTEMBER  1994 


1.  J.O.  Jevtic  and  R.  Lee,  “Techniques  for  Improving  the  Accuracy  of  the  Measured 
Equation  of  Invariance,”  in  preparation  to  be  submitted  to  Journal  of  Electromagnetic 
Waves  and  Application. 

2.  R.G.  Rojas,  “Generalized  Impedance/Resistive  Boundary  Conditions  for  a  Planar  Ho¬ 
mogeneous  Chiral  Slab.” 

3.  M.F.  Otero  and  R.G.  Rojas,  “Two  Dimensional  Green’s  Function  for  an  Impedance 
Wedge.” 

4.  M.F.  Otero  and  R.G.  Rojas,  “Scattering  of  a  Complex  Material  Body  of  Arbitrary 
Shape  in  the  Presence  of  an  Impedance  Wedge.” 

5.  Y.M.  Hwang,  R.G.  Kouyoumjian  and  M.  Hsu,  “Asymptotic  Analysis  of  the  Scattering 
by  a  Wedge  Illuminated  by  a  Field  with  a  Rapid  Spatial  Variation.’ 

6.  R.J.  Burkholder  and  P.H.  Pathak,  “A  Generalized  Ray  Expansion  for  Computing  the 
EM  Fields  Radiated  by  an  Antenna  in  a  Complex  Environment.” 

7.  M.  Hsu,  P.H.  Pathak  and  C.W.  Chuang,  “Analysis  of  the  Asymptotic  HF  EM  Coupling 
Between  Sources  Anywhere  in  the  Vicinity  of  a  Circular  Cylinder.” 

8.  P.H.  Pathak,  R.J.  Burkholder  and  P.  Rousseau,  “On  the  Question  of  Causality  Associ¬ 
ated  with  the  Inversion  into  Time  Domain  of  Ray  Fields  that  Pass  Through  Caustics.” 

9.  G.  Zogbi,  R.J.  Burkholder  and  P.H.  Pathak,  “An  Efficient  Planar  Antenna  Near  and 
Far  Field  Analysis  using  Gaussian  Aperture  Elements.” 

10.  P.H.  Pathak,  A.  Nagamune  and  R.G.  Kouyoumjian,  “An  Analysis  of  Compact  Range 
Measurements.” 


5 


JSEP  RELATED  CONFERENCES/ORAL  PRESENTATIONS 
SEPTEMBER  1993  TO  SEPTEMBER  1994 

1.  Y.S.  Choi-Grogan,  R.  Lee,  K.  Eswar  and  P.  Sadayappan,  “The  Performance  of  a  Parti¬ 
tioning  Finite  Element  Method  on  the  Touchstone  Delta,”  accepted  to  the  10th  Annual 
Review  of  Progress  in  Applied  Computational  Electromagnetics,  Monterey,  CA,  March 
1994. 

2.  A.  Bataineh,  R.  Lee  and  F.  Ozguner,  “Electrical  Characterization  of  High-Speed  In¬ 
terconnects  with  a  Parallel  Three-Dimensional  Finite  Difference  Time-Domain  Algo¬ 
rithm,”  High  Performance  Computing  ’94,  Lojolla,  California,  April  1994. 

3.  M.F.  Otero  and  R.G.  Rojas,  “Radiation/Scattering  of  Complex  Material  Body  of  Arbi¬ 
trary  Shape  in  the  Presence  of  a  Material  Loaded  Impedance  Wedge,”  IEEE  APS/URSI 
International  Symposium,  Seattle,  Washington,  June  1994. 

4.  R.G.  Kouyoumjian,  “Roger  and  Me,  ”  IEEE  APS/URSI  International  Symposium, 
Seattle,  Washington,  June  1994. 

5.  P.R.  Rousseau  and  P.H.  Pathak,  “Time  Domain  Version  of  the  UTD  for  a  Curved 
Wedge,”  IEEE  APS/URSI  International  Symposium,  Seattle,  Washington,  June  1994. 

6.  C.H.  Chuang,  P.H.  Pathak  and  R.J.  Burkholder,  “Hybrid  Analysis  of  EM  Scattering 
by  Large  Open  Cavities  including  Effects  of  Aperture  Shape  and  Interior  Obstacle,” 
IEEE  APS/URSI  International  Symposium,  Seattle,  Washington,  June  1994. 

7.  E.H.  Newman  and  I.  Tekin,  “An  Overview  of  the  Method  of  Moments  Analysis  of 
Large  Systems,”  Board  of  Mathematical  Sciences  National  Research  Council  Sympo¬ 
sium  on  Large-Scale  Structures  in  Acoustics  and  Electromagnetics,  National  Academy 
of  Sciences,  Washington,  DC,  September  26-27,  1994. 

8.  M.  Hsu  and  P.H.  Pathak,  “Hybrid  Analysis  (MM-UTD)  of  EM  Scattering  by  Finned 
Objects,”  IEEE  APS/URSI  International  Symposium,  Seattle,  Washington,  June  1994. 

9.  R.J.  Burkholder,  “A  Progressive  Physical  Optics  Algorithm  for  Computing  the  EM 
Scattering  by  Jet  Inlet  Cavities,”  IEEE  APS/URSI  International  Symposium,  Seattle, 
Washington,  June  1994. 

10.  G.  Manara,  P.  Nepa  and  R.G.  Kouyoumjian,  “Three  Dimensional  Solution  for  the 
Diffraction  of  an  Inhomogeneous  Plane  Wave  by  a  Wedge,”  IEEE  APS/URSI  Interna¬ 
tional  Symposium,  Seattle,  Washington,  June  1994. 

11.  M.  Hsu,  P.H.  Pathak  and  H.  Tseng,  “Hybrid  (MM-UTD)  Analysis  of  EM  Scattering 
by  Large  Convex  Objects  with  Appendages,”  Applied  Computational  Electromagnetics 
Society  (ACES)  Conference,  Monterey,  CAj  March  1994. 
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12.  R.J.  Burkholder,  P.R.  Rousseau  and  P.H.  Pathak,  “  A  Hybrid  Approach  for  Computing 
the  EM  Scattering  from  Complex  Terminations  Inside  Large  Open  Cavities,”  Applied 
Computational  Electromagnetics  Society  (ACES)  Conference,  Monterey,  CA,  March 
1994. 

13.  G.  Zogbi,  H.T.  Chou  and  P.H.  Pathak,  “Reflection  and  Transmission  of  Well  Fo¬ 
cussed  General  Astigmatic  EM  Gaussian  Beams,”  International  Symposium  on  Electro¬ 
magnetic  Environments  and  Consequences,  EURO  EM  (NEM/HPEM)’94,  Bordeaux, 
Prance,  May-June  1994.  (Invited) 

14.  P.R.  Rousseau  and  P.H.  Pathak,  “Development  of  TD-UTD  and  Its  Modiflcations  for 
Analyzing  the  Transient  Scattering  from  Curved  Wedge  Configurations,”  5th  Sympo¬ 
sium  on  Mathematical  Methods  in  Electromagnetic  Theory  (MMET*94),  Kharkov, 
Ukraine,  September  1994.  (Invited) 

15.  J.O.  Jevtic  and  R.  Lee,  “On  the  Choice  of  Metrons  in  the  MEI  Method,”  2nd  Workshop 
on  Finite  Element  Methods  in  Electromagnetic  Wave  Problems,  Siena,  Italy,  May  1994. 

16.  R.  Lee  and  J.O.  Jevtic,  “Several  Approaches  to  Improving  the  Solution  in  the  Measured 
Equation  of  Invariance,”  Joint  AP-S  Symposium  and  URSI  Radio  Science  Meeting, 
Seattle,  WA,  June  1994.  (Invited) 

17.  T.T.  Chia,  R.J.  Burkholder  and  R.  Lee,  “Implementing  the  Hybrid  Ray-FDTD  Method 
for  Computing  the  RCS  of  3-D  Open-Ended  Waveguide  Cavities,”  Joint  AP-S  Sympo¬ 
sium  and  URSI  Radio  Science  Meeting,  Seattle,  WA,  June  1994. 

18.  U.  Pekel  and  R.  Lee,  “A  Three-Dimensional  Bymoment  Method  Approach  for  the 
Analysis  of  Electromagnetic  Scattering,”  Joint  AP-S  Symposium  and  URSI  Radio 
Science  Meeting,  Seattle,  WA,  June  1994. 
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JSEP  RELATED  AWARDS,  PH.D.  DISSERTATIONS  AND  M.SC.  THESES 
SEPTEMBER  1003  TO  SEPTEMBER  1004 

Awards: 

1.  L.M.  Chou,  1993  ESL  Best  Dissertation  Award,  “A  Novel  Hybrid  Full- Wave  Analysis 
Method  for  Planar  Transmission  Lines  Embedded  in  Multilayered  Dielectrics  —  The 
WH/GSMT,”  Supervisor:  R.G.  Rojas. 

2.  Prof.  Pathak  has  received  the  1994  Lumley  Research  Award  from  the  College  of  Engi¬ 
neering  at  The  Ohio  State  University.  This  award  is  based  on  the  previous  four  years 
of  work. 

Dissertations: 

1.  George  Zogbi,  “Reflection  and  Diffraction  of  General  Astigmatic  Gaussian  Beams  from 
Curved  Surfaces  and  Edges,  ”  June  1994. 

2.  Tse-Tong  Chia,  “Application  of  FDTD  in  Hybrid  Methods  for  EM  Scattering  Analysis 
of  Cavities  with  Complex  Terminations,”  June  1994. 

3.  M.E.  Peters,  “Method  of  Moments  Analysis  of  Anisotropic  Artificial  Media  Composed 
of  Dielectric  Wire  Objects,”  March  1994. 

Theses: 

1.  Ty  Barkdoll,  “Finite  Element  Analysis  of  Bodies  of  Revolution  using  MEI  Boundary 
Conditions,”  September  1993. 

2.  Yung  Shirley  Choi-Grogan,  “A  Sequential  and  Parallel  Implementation  of  a  Partition¬ 
ing  Finite  Element  Technique  for  Electromagnetic  Problems,”  December  1993. 

3.  Jovan  Jevtic,  “An  Analysis  of  the  Measured  Equation  of  Invariance,”  May  1994. 
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Radio  Science,  Volume  29,  Number  2,  Pages  465-481,  March-April  1994 


Efficient  numerical  and  closed-form  asymptotic  representations 
of  the  dyadic  aperture  Green’s  function  for  material-coated 
ground  planes 

Gary  A.  Somers 

Lincoln  Laboratory,  Massachusetts  Institute  of  Technology,  Lexington 

Prabhakar  H.  Pathak 

The  Ohio  State  University,  ElectroScience  Laboratory,  Columbus 


Abstract.  This  paper  provides  both  numerical  and  closed-form  uniform  asymptotic 
representations  of  the  aperture  Green’s  function  which  determines  the  magnetic  field  • 
due  to  a  point  magnetic  current  source  where  both  the  source  and  field  points  lie  on 
the  ground  plane  which  is  covered  by  a  planar  material  slab.  The  closed-form 
asymptotic  representation  which  explicitly  contains  effects  of  space,  surface,  and  leaky 
waves  is  shown  to  agree  with  the  numerical  evaluation  of  the  exact  integral 
representation  for  source  and  field  point  separations  which  are  only  a  few  tenths  of  a 
free  space  wavelength.  An  extended  envelope  extraction  technique  allows  an  efficient 
evaluation  of  the  exact  integral  representation  for  source  and  observation  separations 
too  small  for  the  asymptotic  solution  to  be  valid.  The  asymptotic  solution  is  uniform  in 
that  it  remains  valid  across  all  of  the  relevant  surface  and  leaky  wave  transition 
regions.  The  computation  time  for  this  asymptotic  result  is  2-3  orders  of  magnitude 
faster  than  the  numerical  evaluation  of  the  exact  integral  form,  and  additionally  its 
computation  time  is  independent  of  the  source  and  field  point  separations.  This 
aperture  Green’s  function  is  very  useful  when  computing  the  mutual  coupling  elements 
of  the  method  of  moments  admittance  matrix  for  an  array  of  slots  in  a  material  covered 
ground  plane. 


1.  INTRODUCTION 

There  is  increasing  interest  in  analyzing  electri¬ 
cally  large  antennas  that  are  embedded  in  a  material 
coating.  A  large  array  of  slots  in  a  material-coated 
ground  plane  fall  into  this  category.  When  using  an 
integral  equation  formulation  for  the  surface  field, 
which  is  solved  via  the  method  of  moments  (MM), 
it  is  necessary  to  calculate  the  mutual  coupling 
between  all  of  the  slots  in  the  array.  These  mutual 
coupling  calculations  typically  involve  Sommerfeld 
integrals  whose  brute  force  numerical  evaluation 
can  be  a  formidable  task  especially  when  the  slots 
are  separated  by  large  electrical  distances.  In  this 
paper,  two  very  efficient  representations  of  the 
dyadic  aperture  Green’s  function  are  presented 
which  facilitates  these  mutual  admittance  calcula¬ 
tions  that  otherwise  might  have  been  intractable  for 
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large  arrays.  The  numerical  solution  implemented 
here  complements  the  closed  form  uniform  asymp¬ 
totic  solution  in  that  it  is  accurate  and  efficient  in  the 
small  region  surrounding  the  source  point  where  the 
asymptotic  solution  fails.  Additionally,  the  numer¬ 
ical  solution  serves  as  a  reference  solution  for  the 
asymptotic  formulation.  These  solutions  are  useful 
for  analyzing  slots  in  a  ground  plane  which  are 
excited  by  striplines  (or  other  waveguiding  struc¬ 
tures)  in  multilayered  substrates  below  the  ground 
plane. 

The  aperture  dyadic  Green’s  function  presented 
here  provides  the  magnetic  fields  on  a  ground  plane 
with  a  material  coating  caused  by  a  unit  magnetic 
current  point  source  which  is  also  on  the  ground 
plane  as  shown  in  Figure  1.  It  is  assumed  that  the 
ground  plane  is  a  perfect  electric  conductor  (PEC), 
and  the  material  coating  is  isotropic  and  homoge¬ 
neous.  It  is  also  assumed  that  the  material  above  the 
slab  is  free  space.  Recently,  there  has  been  activity 
in  the  related  area  of  calculating  the  electric  current/ 
electric  field  Green’s  function,  where  both  the 
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Figure  l._  Material-coated  ground  plane  with  magnetic 
current,  M,  at  f  and  observation  point  at  r. 


source  and  observation  points  are  on  the  material- 
air  interface  of  the  grounded  planar  material  slab 
[Barkeshli  et  al.,  1990;  Marin  et  al.,  1989;  Marin 
and  Pathak,  1992].  All  of  these  papers  share  a 
common  formulation  in  that  they  calculate  poten¬ 
tials  which  then  need  to  be  diflFerentiated  either 
numerically  or  analytically  to  find  the  fields.  In  this 
work  the  Green’s  function  is  determined  directly 
which  eliminates  the  need  to  differentiate  an  asymp¬ 
totic  series.  In  general,  it  is  preferable  to  develop  an 
asymptotic  approximation  of  a  differentiated  ex¬ 
pression  rather  than  differentiating  an  asymptotic 
approximation.  Additionally,  in  the  previous  re¬ 
lated  approach  [Barkeshli  et  al.,  1990;  Marin  et  al., 
1989;  Marin  and  Pathak,  1992]  the  transverse  elec¬ 
tric  (TE)  and  transverse  magnetic  (TM)  potentials 
were  coupled;  in  contrast,  this  work  keeps  them 
uncoupled.  The  coupled  solution  leads  to  more 
complicated  integrand  topologies  in  the  spectral 
plane  on  which  the  asymptotic  approximations  are 
performed  when  compared  to  the  uncoupled  topol¬ 
ogies.  It  is  preferable  to  deal  with  the  simplest 
possible  topology  from  an  accuracy  perspective 
when  performing  an  asymptotic  analysis.  One  final 
major  difference  between  the  analysis  presented 
here  and  by  Barkeshli  et  al.  [1990],  Marin  et  al. 
[1989],  and  Marin  and  Pathak  [1992]  is  that  in  the 
present  asymptotic  approximation,  the  influence  of 
any  poles  on  the  saddle  point  contribution  is  com¬ 
pensated  by  using  the  modified  Pauli-Clemmow 


technique  [Clemmow,  1950]  which  leads  to  a  multi¬ 
plicative  correction  factor  (transition  function), 
whereas  [sec  Barkeshli  et  al.,  1990;  Marin  et  al., 
1989;  Mari,"  and  Pathak,  1992]  the  Van  der  Waer- 
•  den  Method  [Van  der  Waerden,  1951;  Felsen  and 
Marcuvitz,  1973]  was  applied  which  results  in  an 
additive  correction  term.  The  multiplicative  correc¬ 
tion  is  preferable  because  calculations  are  not  very 
sensitive  to  the  computational  precision.  When  the 
correction  factor  is  additive  and  a  surface  or  leaky 
wave  pole  (LWP)  is  close  to  the  saddle  point,  then 
the  solution  requires  one  to  take  the  difference  of 
two  very  large  numbers  which  numerically  might 
not  produce  correct  results.  The  advantage  to  using 
the  Van  der  Waerden  method  is  that  if  there  are 
many  poles  that  need  to  be  accounted  for,  it  is  a 
natural  extension  to  simply  compensate  for  each  of 
the  poles.  When  using  the  modified  Pauli-Clemmow 
technique,  a  new  correction  factor  must  be  derived 
for  each  multiplicity  of  poles.  It  is  important  to  note 
that  in  practice,  it  is  unlikely  that  it  will  be  neces¬ 
sary  to  compensate  for  more  than  three  poles, 
although  this  paper  provides  the  expressions  for  an 
arbitrary  number  of  poles.  It  will  be  shown  that 
these  expressions  are  easily  calculatable,  highly 
accurate,  and  efficient. 

This  paper  is  organized  as  follows.  Section  2 
presents  the  formulation  of  dyadic  Green’s  function 
in  terms  of  Sommerfeld  integrals.  Section  3  pre¬ 
sents  a  scheme  referred  to  as  the  envelope  extrac¬ 
tion  technique  to  numerically  evaluate  the  exact 
solution  in  terms  of  the  four  Sommerfeld  integrals. 
It  is  noted  once  again  that  the  numerical  solution  is 
used  not  only  as  a  reference  to  validate  the  accu¬ 
racy  of  the  asymptotic  solution,  but  also  to  comple¬ 
ment  the  asymptotic  solution  in  the  small  region 
where  the  latter  becomes  invalid.  The  uniform 
asymptotic  evaluation  of  the  Sommerfeld  integrals 
will  be  presented  in  section  4  where  some  common 
cases  will  be  addressed  along  with  the  general 
formulation.  Finally,  section  5  will  show  some 
numerical  results  displaying  the  accuracy  of  the 
asymptotic  solution,  the  computation  times  and 
some  concluding  remarks.  An  e'"'  time  dependence 
is  assumed  and  suppressed  in  the  following  analysis. 

2.  Dyadic  Green’s  Function  Formulation 

If  we  consider  an  arbitrarily  directed  impressed 
magnetic  surface  current  density,  M^,  on  the 
ground  plane  at  any  point  f*  over  some  region  s,  and 
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we  wish  to  determine  the  magnetic  field  H  on  the 
ground  plane  at  any  point  F,  the  following  expres¬ 
sion  can  be  used  [Tai,  1973]: 


H(r)  = 


G(r-,  n  •  Ms(n  ds'. 


(1) 


The  material  coating  on  the  PEC  and  free  space 
above  it  have  the  constitutive  parameters  ej,  /H)  and 
eq,  /liq,  respectively.  Since  the  normal  magnetic 
field  vanishes  on  a  PEC,  then  by  reciprocity,  the 
normal  magnetic  current  does  not  radiate.  We  can 
express  the  Green’s  function  in  a  rectangular  (x,  y, 
z)  basis  as  follows: 


G{r-,  n  = 


-G,Ar,n  G,y(r-,n  0-| 

Gy^(r-  F)  Gyy(r;  F)  0 
0  0  0 


(2) 


makes  with  the  positive  x  axis.  These  Sommerfeld 
integral  terms  are  only  a  function  of  one  variable  p; 
therefore  they  can  be  calculated  once  and  stored  in 
a  “look-up  table”  and  later  recalled  to  construct  the 
Green’s  functions.  This  is  particularly  useful  when 
analyzing  a  large  two-dimensional  array. 

The  Sommerfeld  integrals  are  given  by 

TM  r*  tDE  1  A(z) 

and 

TE  r“ 


where 


The  first  index  refers  to  the  polarization  of  the 
magnetic  field  component,  and  the  second  index 
refers  to  the  direction  of  the  magnetic  current 
component  in  rectangular  coordinates.  Each  of  the 
four  nonzero  elements  of  the  dyadic  Green’s  func¬ 
tion  above  can  be  expressed  in  terms  of  TE^  and 
TMj  Sommerfeld  integrals: 

G,Ar;r')  =  ^U™(p)-/l^(p)) 

cos  (20) 

+  (/2™(P)  +  /p(p))  (3) 

sin  (20) 

G^.(F;  F' )  =  GyAF;  F' )  =  (/J%)  +  /p(p))  (4) 


,(f;r-')  =  -(/™(p)-/TE(p)) 


cos  (20) 


,TM 


where 


'(p)  +  /p(p)) 

(5) 

p  =  |r--r-'| 

(6) 

0  =  tan  '( -  . 

\(r-F')-x 

(7) 

In  other  words,  p  is  the  distance  between  the  source 
and  observation  points,  and  0is  the  angle  the  vector 
from  the  source  point  to  the  observation  point 


A(z)  =  -eq/:!''  cos  -  z)) 

sin  (k['\h-z)) 

(10) 

B(z)  =  M,Af  cos  ()t''*(/;  -  z)) 

+  sin  (k'^'\h  -  z)) 

(11) 

TM  =  eAz’  cos  (^y’/i)  +jEoi(*”  sin  {kl% 

(12) 

TE  =  cos  (kl'^h)  +jp,]k^^^  sin  {k[^^h), 

(13) 

where  z  is  the  height  off  the  ground  plane  of  the 
observation  point,  ft  is  the  material  thickness, 
and  are  the  z-directed  propagation  constants  in 
free  space  and  in  the  material,  respectively.  Since 
corresponds  to  a  semiinfinite  region  (z  =  /i  -» 
»),  it  contains  a  branch  cut  and  since  corre¬ 
sponds  to  a  finite  region  (z  =  0  -»  h),  it  is  single 
valued.  These  properties  can  be  verified  by  the 
integrands  of  (8)  and  (9). 

=  (*o^-^2)1/2;  (14) 

and 


=  (15) 

where  kg  and  are  the  propagation  constants  in  free 
space  and  in  the  material,  respectively.  For  the  loss¬ 
less  case,  there  will  be  at  least  one  surface  wave  pole 
(SWP)  on  the  real  axis  of  the  Sommerfeld  plane  (zeros 
of  TM  and/or  TE).  To  satisfy  the  radiation  condition, 
the  contour  must  be  skirted  above  the  poles. 
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As  indicated  previously,  there  will  be  two  meth¬ 
ods  presented  to  accurately  and  efficiently  compute 
the  Sommerfeld  integrals  in  (8)  and  (9).  One  method 
is  a  numerical  method  which  works  well  for  very 
small  p  and  hence  complements  the  asymptotic 
solution  since  the  alsymptotic  solution  is  accurate 
for  relatively  small  to  any  large  value  of  p. 


3.  Numerical  Evaluation  of  the  Exact 
Solution  for  the  Green’s  Dyacic 
(for  p  Small) 

In  this  section  a  scheme  to  enhance  the  conver¬ 
gence  of  the  Sommerfeld  integrals  which  comprise 
the  Green’s  dyadic,  G,  as  indicated  in  (3)-(5)  will  be 
presented.  This  numerical  solution  complements 
the  closed-form  uniform  asymptotic  solution  (devel¬ 
oped  in  the  next  section)  in  a  region  for  p  small, 
which  is  where  the  asymptotic  solution  fails.  It  is 
noted  that  the  asymptotic  solution  is  useful  as  long 
as  the  separation  between  the  source  and  field 
points  in  greater  than  only  a  few  tenths  of  a  free 
space  wavelength.  This  numerical  solution  is  very 
similar  to  what  was  done  by  Jackson  and  Alexopou- 
los  [1986]  with  the  exception  that  for  the  TE  type 
integrals  in  this  paper,  two  terms  were  subtracted 
from  the  integrand  which  leads  to  faster  convergence 
than  if  one  term  were  subtracted  as  done  previously. 
A  related  technique  was  used  by  Pozar  [  1984, 1985]  to 
invert  the  spectral  integral  of  the  MM  mutual  imped¬ 
ance,  whereas  in  this  paper  the  spectral  integrals  are 
inverted  to  produce  the  spatial  domain  representation 
of  the  various  terms  of  the  Green’s  function. 

The  behavior  of  the  integrands  for  large  ^  in  (8) 
and  (9)  are  given  below: 

integrand  of  ~  |-i- 


•  cos  ~  i  (17) 

If  these  integrands  are  examined  for  the  case  when 
the  observation  point  is  on  the  ground  plane,  that  is, 
when  z  =  0,  (as  is  necessary  in  the  mutual  coupling 


problem  for  slots),  then  the  integrands  no  longer 
exhibit  an  exponential  decay.  The  TM  integrands 
have  an  algebraic  decay,  and  the  TE  inte¬ 
grands  show  an  algebraic  growth.  The  ip  term 
in  the  argument  of  the  cosine  functions  leads  to 


Normolized  Wavenumber 


Normalized  Wovenumber 

Figure  2.  Integrand  of  /™  for  (a)  p  =  0.2Ao,  (b)  p  = 
0.5Ao  and  (c)  p  =  2.0Ao.  The  material  slab  is  characterized 
by  Sr  =  2.2,  p,  =  1.0,  and  a  height  of  0.2Ao. 
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faster  integrand  oscillation  when  the  separation  be¬ 
tween  the  source  and  observation  points  p  becomes 
large.  These  TM  and  TE  integrands  are  shown  in 
Figures  2  and  3,  respectively,  for  various  values  of  p. 
These  integrands  have  such  poor  convergence  (and 
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Normolized  Wavenumber 

Figure  3.  Integrand  of  /p  for  (a)  p  =  0.2Ao,  (b)  p  = 
O.SAo  and  (c)  p  =  2.0Ao.  The  materia!  slab  is  characterized 
by  E,  -  2.2,  p.,  =  I.O,  and  a  height  of  0.2Ao. 


divergence  in  the  TE  case  if  z  =  0)  properties  because 
the  Sommerfeld  representation  is  a  z-propagating  rep¬ 
resentation,  and  for  this  case  there  is  no  z  separation 
between  the  source  and  observation  points.  If  a  brute 
force  numerical  integration  was  to  be  performed,  then 
it  would  be  wise  to  transform  the  integrand  to  a 
p-propagating  radial  representation  [Felsen,  1973]. 
However,  since  this  solution  is  also  intended  to  be 
used  for  small  p,  a  radial  representation  would  not  be 
a  wise  choice.  The  envelope  extraction  technique  is 
applied  here  [Jackson  and  Alexopoulos,  1986].  In  this 
technique  a  large  argument  form  of  the  integrand  is 
found  that  can  be  integrated  analytically.  This  result  is 
then  added  to  one  in  which  the  large  argument  form  is 
subtracted  from  the  original  integrand  leaving  a 
quickly  converging  difference  to  be  integrated  numer¬ 
ically.  The  zeros  of  TM  and  TE  are  spectral  poles 
which  represent  surface  wave  contributions.  In  gen¬ 
eral,  there  can  be  a  multitude  of  zeros  of  both  TM  and 
TE.  Analogous  to  the  envelope  extraction  technique, 
the  singularities  can  be  subtracted  out  of  the  integrand 
and  replaced  analytically  since  the  pole  contribution 
can  be  integrated  in  closed  form.  It  is  only  necessary 
to  subtract  the  singularities  in  the  region  surrounding 
them,  say  from  to  ^2  such  that  <  ^2  for 

all  pole  locations  n  =  1,  2,  •  •  • 


CUE)  A(z  -  0) 
Tp"  TM 


+ 


WE  I  A{z  =  0) 
TM 


Res  {FTMif™)}] 

^  -  C  1 


d^ 


f“  [wEl 

A(z  =  0) 

L 

TM 

JWB] 


+  2  Res  (EtmO) 


In 


/:TM  _  , 


■^J 


(18) 


where  n  is  summed  over  the  finite  number  of  TM 
SWPs  that  are  solutions  to  TM(^™)  =  0.  The 
residue  is  given  by 
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Res{FTM(0>^ 


v<oe\  A{z  =  0) 


a 

^  —  TM 

H 


(19) 


Figure  4  shows  an  example  of  the  modified  inte¬ 
grand  for  the  cases  shown  in  Figure  2,  with  the 
envelope  and  singularity  extraction  techniques  ap¬ 
plied.  Note  the  scale  change  of  both  axes.  These 
new  integrands  are  less  oscillatory  and  decay  fast; 
hence  they  are  much  easier  to  numerically  inte¬ 
grate.  There  is  only  a  single  SWP  for  this  configu¬ 
ration.  The  spike  in  the  plots  shows  that  there  was 
not  a  perfect  cancellation  at  the  singularity.  Numer¬ 
ically,  this  is  not  a  problem  because  a  numerical 
integration  technique  such  as  Gauss  quadrature  that 
does  not  integrate  at  the  endpoint  of  an  interval  can 
be  used  to  guarantee  that  the  pole  location  will  not 
be  chosen  for  evaluation. 

The  TE  Sommerfeld  integrals  will  be  handled  in 
much  the  same  manner  as  the  TM  Sommerfeld 
integrals.  Mathematically,  in  the  present  form  the 
TE  integrands  are  divergent  for  2  =  0  and  are 
therefore  meaningless.  This  does  not  mean  that  the 
proper  information  is  not  contained  in  these  inte¬ 
grals;  rather,  it  is  indicative  of  the  failure  of  the 
Sommerfeld  representation  for  z  =  0.  If  we  consider 
the  fields  at  a  small  distance  above  the  ground 
plane,  perform  the  envelope  subtraction,  and  take 
the  limit  as  this  distance  goes  to  zero,  useful  expres¬ 
sions  are  obtained.  A  minor  difference  between  the 
TE  and  TM  integrals  is  that  it  is  possible,  and 
beneficial,  to  subtract  2  orders  of  the  large  argument 
form  of  the  integrand  for  the  TE  case: 


B(z  =  0) 
iofij  TE  ^ 


/: 


A:'”  B(z  =  0) 


O 

J 

Wfli 


(oHi  TE 


Jk} 


2o)ni 


J  {3]  jk]  1 

- -^^^-1- - 

£|>M1  P  2w/xi  p 

+  X  Res  {FTE(f Jf)} 


-  Q 


- 


(20) 


where  n  is  summed  over  the  finite  number  (if  any)  of 
TE  SWPs,  which  are  solutions  to  TE(^J^^)  =  0, 
and 


Res  {FTE(f  J.^)}  = 


B(z  =  0) 


W/JL]  d 


TE 


■  (21) 


It  is  a  common  occurrence  for  substrates  that  are 
electrically  thin  not  to  have  any  TE  surface  wave 
poles.  For  this  case,  (20)  takes  the  following  form: 


Biz  =  0) 
TE 


e 


jk] 

2(i>p] 


J[^i^p)d^- 


J 


Jk]  1 

+ - 

2a>p]  p 


(22) 


The  integrand  for  the  latter  case  is  shown  in  Figure 
5,  indicating  that  it  is  slowly  oscillating  and  quickly 
convergent;  thus  it  is  in  a  form  amenable  to  numer¬ 
ical  integration. 
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Res  {FTE(fr)}l 


TE 


d^ 


4.  Uniform  Asymptotic  Closed  Form 
Evaluation  of  the  Green’s  Dyadic  (for  p 
Relatively  Small  to  Large) 

In  this  section,  closed-form  uniform  asymptotic 
solutions  for  the  four  Sommerfeld  integrals  in  (8)- 
(9)  are  presented.  This  asymptotic  solution  comple¬ 
ments  the  numerical  solution  of  the  previous  sec¬ 
tion  since  the  asymptotic  solution  is  accurate  and 
efficient  for  p  greater  than  a  few  tenths  of  a  free 
space  wavelength  below  which  the  numerical  solu- 


SOMERS  AND  PATHAK:  GROUNDED  SLAB  APERTURE  GREEN’S  DYAD 


471 


tion  is  useful.  These  results  can  be  substituted  in 
(3)-(5)  to  directly  arrive  at  the  closed-form  asymp¬ 
totic  approximation  for  the  Green’s  dyadic  G  of  (2). 
The  incentive  for  using  this  asymptotic  solution  is 


Figure  4.  Integrand  of  with  the  envelope  and  singu¬ 
larity  extracted  as  indicated  in  (18)  for  (a)  p  =  0.2Ao,  (b) 
p  =  O.SAfl  and  (c)  p  =  2.0Ao.  The  parameters  e,  =  2.2,  /a, 
=  1.0,  and  h  =  0.2Ao. 
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Figure  5.  Integrand  of  with  the  envelope  extraction 
as  indicated  in  (22)  for  (a)  p  =  0.2Ao,  (b)  p  =  O.SAq  and  (c) 
p  =  2.OA0.  The  parameters  e,  =  2.2,  10.  and  h  = 

0.2Ao. 


that  it  is  extremely  computationally  efficient  and 
typically  remains  accurate  even  for  source  and 
observation  point  separations  of  only  a  few  tenths 
of  a  free  space  wavelength  because  of  its  uniform 
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nature.  All  of  the  relevant  transition  effects  associ¬ 
ated  with  the  confluence  of  surface  and  leaky  do¬ 
mains  with  the  space  waves  directly  radiated  by  the 
source  are  incorporated  into  the  solution. 

Equations  (8)  and  (9)  are  expressed  in  terms  of  a 
semiinfinite  contour  which  is  not  convenient  for 
function  theoretic  manipulations  of  integrals.  A 
more  appropriate  form  is  given  below; 

TM  1  r*  (2) 

"2  J  jO|(fp)  f  (23) 

TE  If®  (21 

/j■}(p)  =  2J  fTEW|0|(fp)  f  (24) 

where 

we  I  A(z  =  0) 


^TE  = 


Biz  =  0) 
w/jL]  TE 


(26) 


Note  that  H^\^p)  has  a  branch  cut  on  the  negative 
real  |  axis,  and  has  a  double  pole  at  the 

origin. 

The  modified  Pauli-Clemmow  technique  [Clem- 
mow,  1950],  which  is  used  here,  corrects  for  the 
influence  of  the  SWPs  and  LWPs  on  the  saddle 
point  by  multiplying  the  conventional  space  wave 
solution  by  the  appropriate  transition  function.  The 
transition  function  that  should  be  used  depends  on 
the  number  of  poles  that  are  in  the  saddle  point 
region.  For  most  problems  of  practical  interest,  it 
appears  sufficient  to  consider  the  influence  of  one  to 
three  poles.  The  number  of  poles  that  need  to  be 
considered  is  discussed  by  Marin  et  al.  [1990]. 

Unlike  the  numerical  method  for  evaluating  G 
discussed  in  section  2,  where  it  is  only  necessary  to 
be  concerned  with  the  proper  (or  top)  Riemann 
sheet  of  the  Sommerfeld  plane,  the  asymptotic 
method  implicitly  incorporates  the  features  of  the 
improper  (or  bottom)  Riemann  sheet,  and  thus  one 
must  be  cognizant  of  its  topology.  It  is  visually 
instructive  (although  not  necessary)  to  map  both  the 
improper  and  proper  Riemann  sheets  of  the  = 
{kl  -  ^^)'''^  function  onto  the  angular  spectral,  y, 
plane  by  the  following  mapping: 


Figure  6.  Angular  spectral  plane  with  the  real  axis  of 
the  Sommerfeld  plane  mapped  into  F,  and  the  steepest 
descent  path  is  indicated  by  Fsdp’  The  saddle  point  is  y^. 


•W|0|(/:oP  sin  y)  sin  y  cos  y  dy 


(28) 


^0  f  (2) 

(p)  =  Y  j  Fje  //|Oj(A:oP  sin  y)  sin  y  cos  y  dy 


(29) 


where  F  is  the  y-plane  image  of  ^plane  real  axis  on 
the  proper  sheet  under  the  mapping  in  (27)  (see 
Figure  6).  The  integrals  in  (28)  and  (29)  are  exact 
representations  of  the  solution.  To  facilitate  the 
asymptotic  approximation,  the  Hankel  functions 
are  replaced  by  their  two-term  large  argument  form. 
This  is  justified  (with  hindsight)  because  the  saddle 
point  is  sufficiently  far  from  the  origin: 


(30) 


I  =  *0  sin  y  (27) 

The  integrals  in  (23)  and  (24)  under  the  above 
mapping  take  the  following  form; 


In  light  of  the  asymptotic  approximations  for  the 
Hankel  functions  above,  it  is  clear  that  the  inte¬ 
grands  of  (28)  and  (29)  possess  a  saddle  point,  y^,  at 
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■y^  =  77/2.  The  original  contour,  T,  can  be  deformed 
to  the  steepest  descent  path  (SDP)  shown  in  Figure 
6.  All  of  the  surface  wave  poles  are  captured  in  the 
contour  deformation.  The  relevant  poles  all  lie  in 
the  region  of  the  saddle  point  along  the  line  Re  (-y)  = 
ttH  (for  the  lossless  case).  The  poles  in  this  region 
that  lie  in  the  upper  half  of  the  -y  plane  (Im  {y}  >  0), 
are  surface  wave  poles  and  those  which  lie  in  the 
lower  half  of  the  •)^plane  (Im  {y}  <  0),  are  leaky 
wave  poles.  The  surface  wave  poles  were  mapped 
from  the  proper  sheet  of  the  Sommerfeld  |  plane, 
and  the  leaky  wave  poles  were  mapped  from  the 
improper  sheet  of  the  Sommerfeld  ^  plane.  The 
Sommerfeld  integrals  can  be  expressed  as  follows: 

.™  ,  ,  ,  .  V'  „  ™  snp,™ 

/ji|(p) - 2ttj  ^  ^ 


/|||(P)  ~  -iTTj  2  Resjfj  +  (32) 


where 


TM  f  +  1  *0  /  2 


sw;  r,  .  r  ^ 


'{!) 


2  X’nkap 


1/2 


■I 


1 

-15 


Ftm  y  cos  y  1 .  , 

rsDP  \  sin  yf 


■  exp  (-Jkgp  sin  y)  dy 


(33) 


sdp^te  f  +  l  ^/_2_ 

1~J  2  [nkop 


1/2 


•^SDP 


Fte  sin'^^  y  cos  yl  1  + 


^  M  \ 

J 

Mqp  sin  yj 

•  exp  (-jkgp  sin  y)  dy. 


(34) 


The  sum  over  n  includes  all  surface  wave  poles.  The 
residues  are  given  by 


„  TM 
Resj.j  - 


""kl  we  I  A(z  =  0) 


2  A'"' 


d 

—  TM 
dy 


(2) 


{kgp  sin  y)  sin  y  cos  y 


(35) 


y  =  yp. 


r.  TE 

ReB{,j  - 


rto  *?’  B(z  -  0) 


2  wfi\ 


p  sin  y)  sin  y  cos  y 


TE 


7=  rpn 


(36) 


Note  that  the  residues  are  extracted  from  the  inte¬ 
gral  before  the  Hankel  function  approximation  of 
(30)  is  implemented;  thus  the  residue  contributions 
are  exact. 

The  standard  form  in  which  to  implement  the 
asymptotic  approximations  of  the  integrals  of  (33)  and 
(34)  can  be  realized  by  mapping  this  y-plane  represen¬ 
tation  to  the  5  plane  by  the  following  transformation: 

sin  y  =  I  -Js^  cos  y  =  ~s(s^  +  2j)'^ 


dy  = 


2j 


(s^  +  2j) 


in 


ds 


(37) 


The  SDP  contribution  of  (31)  and  (32)  can  be 
expressed  in  the  s  plane  by 


SDP 


..TM 
/rn  = 


{^1 


:|l 


2  \ 


in 


nkop 


gJln/4)  g-jk„p 


(i)  exp  {-kops^)  ds 


2  Xyrkop) 


(5)  exp  {-kgps^)  ds 


where 


T{i}(s)  =  -Ftm  (1  s{s^  + 


(38) 


(39) 


(40) 


(41) 
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Since  necessary  to 

perform  a  second-order  asymptotic  evaluation  of 
the  integrals  in  (38)  and  (39).  For  the  case  of  the 
source  and  observation  points  both  on  the  ground 
plane,  which  is  being  addressed  in  this  paper,  the 
SDP  is  fixed  in  the  y  plane  with  the  saddle  point  at 
y  -  7r/2.  When  the  system  parameters  vary,  the 
poles  migrate  and  are  restricted  to  the  line  Re  {y}  = 
7t/2  in  the  vicinity  of  the  saddle  point.  The  only 
point  onto  which  the  poles  can  intersect  the  SDP  is 
the  saddle  point;  thus  it  is  appropriate  to  use  the 
modified  Pauli-Clemmow  technique  [Clemmow, 
1950].  In  practice,  it  is  typically  necessary  to  con¬ 
sider  the  influence  of  one,  two,  or  three  poles 
simultaneously  on  the  SDP  contribution.  These 
cases  are  presented  along  with  the  more  general  N 
pole  case  for  completeness. 

4.1,  Uniform  Asymptotic  Evaiuation  of  an 
Integrai  With  a  Single  Pole  Influence 
on  the  Saddle  Point 

When  implementing  the  modified  Pauli-Clemmow 
technique,  the  function  Tis)  (which  represents  the 

TM  TE 

integrand  in  all  four  cases):  ^nd 

expressed  in  terms  of  a  function  Tgis),  whose  region 
of  analyticity  is  a  disk,  centered  at  5  =  0,  with  a 
radius  equal  to  the  distance  to  the  closest  pole  other 
than  the  pole  being  accounted  for  aX  s  =  Sp]  Sp  can 
be  determined  from  (37): 


T{s)  = 


Ta(s) 

S  Sp 


(42) 


d  wSj^Q  1 

'  d  ko  1 

0) fiQ  cos^  [ik^  -  kl)''-h] 


(45) 


(46) 


and  3'^(x)  is  the  slope  transition  function  [Hwang 
and  Kouyoumjian,  1974]  given  by 

^,(x)  =  2jx[\  -  ^(x)l  (47) 

9'(x)  is  the  wedge  diffraction  transition  function 
[Kouyoumjian  and  Pathak,  1974],  which  is  easily 
calculated  [Pathak  and  Kouyoumjian,  1970]. 

The  argument  of  the  transition  function  can  be 
expressed  as  follows: 


jkop{s™)‘  =  (ko  -  fp™)p 

(48) 

jkopisj^f  =  (ko  -  ^J^)P 

(49) 

where  and  are  the  zeros  of  TM  and  TE, 
respectively.  If  the  pole  under  consideration  is  a 
SWP,  then  the  zero  lies  on  the  top  Riemann  sheet  of 
the  Sommerfeld  plane  and  if  the  pole  corresponds  to 

a  LWP,  then  the  zero  would  lie  on  the  bottom 

rTM| 

Riemann  sheet.  Note  that  since  Re  ^  k^, 

then  the  argument  of  the  transition  function  is  not 
positive,  and  the  following  relationship  is  useful  for 
real  arguments: 


3^(-|;t|)  =  9;*(W)  (50) 


Implementing  the  standard  uniform  asymptotic  re¬ 
duction  of  the  integrals,  one  obtains 


g-jkap 

(AoP)^ 


(43) 


-  {!}^0  Yy  ^sUkoPisJ^?) 


g-Jhp 

8^0  P 

I  (kopf 

where  the  asterisk  indicates  complex  conjugate. 


4.2.  Uniform  Asymptotic  Evaluation 
of  an  Integral  With  Two  Poles 
Near  a  Saddle  Point 

The  development  of  the  case  where  two  poles,  Sp^ 
and  Sp^,  are  close  enough  to  the  saddle  point  to 
influence  its  evaluation  is  similar  to  the  single  pole 
case  developed  in  the  previous  section.  Now  the 
integrand  is  represented  by  the  following  equation: 


T{s)  = 


Ta(s) 

{S  -  5p,)(5  -  Jp,) 


(51) 


Performing  a  partial  fraction  expansion  on  (51)  and 
implementing  an  asymptotic  procedure  similar  to 
the  single  pole  case,  the  uniform  asymptotic  solu¬ 
tion  can  be  found  to  be 


where 
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(53) 


where 

fTM(y.)  =;'  ,  2  '"?2u72  (54) 

(ki  —  kq) 

(^2_^2U/2 

FrB(ys)=j- - —tanlikyk^oy'^h]  (55) 

wpi 

When  evaluating  the  SDP  contribution  that  is  influ¬ 
enced  by  two  or  more  poles,  it  is  necessary  to  use 
the  proper  branch  of  Sp  = 

sp  =  ±e>^^'*^\sp\  (56) 

where  the  plus  sign  is  taken  for  SWPs  and  the  minus 
sign  is  taken  for  LWPs.  The  parameter  Sp  refers  to 


either  a  SWP  or  LWP  that  is  in  the  vicinity  of  the 
saddle  point. 


4.3.  Uniform  Asymptotic  Evaluation 
of  an  Integral  With  Three  Poles 
Near  a  Saddle  Point 

The  three-pole  case  follows  the  same  develop¬ 
ment  as  the  two-pole  case  where  the  integrand  is 
now  represented  by 


Tis)  = 


TAs) 

{s  -  Sp^S  -  Sp^S  -  Sp^ 


(57) 


Performing  the  appropriate  manipulations  and  as¬ 
ymptotic  analyses,  one  obtains  equations  (58)  and 
(59). 
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4.4.  Uniform  Asymptotic  Evaluation 
of  an  Integral  With  N  Poles 
Near  a  Saddle  Point 

If  one  wishes  to  extend  the  asymptotic  analysis 
to  an  electrically  thick  substrate,  then  the  same 
formulation  applies.  In  general,  as  the  material 
coating  becomes  electrically  thicker,  it  becomes 
necessary  to  account  for  more  poles  (SWPs 
and  LWPs)  which  move  into  the  proximity  of 
the  saddle  point.  The  same  technique  is  used 
here  that  was  used  in  sections  4. 1-4.3  to  eval¬ 
uate  an  integral  via  a  uniform  asymptotic  pro¬ 
cedure  which  accounts  for  the  confluence  of  N 
poles  near  a  saddle  point.  The  integrand  is  written 
as  follows: 


N 

n 

m  =  1 

where  N  is  the  number  of  poles  incorporated  into 
the  analysis  of  that  particular  integral.  This  leads  to 
the  closed-form  uniform  asymptotic  approximation 
of  the  SDP  integrals  of  (38)  and  (39): 


5.  Results  and  Conclusions 

In  this  section  we  will  examine  the  accuracy  of 
the  closed-form  uniform  asymptotic  solution  devel¬ 
oped  in  section  4  by  comparing  it  with  the  “exact” 
numerical  solution  which  is  calculated  by  using  the 
envelope  extraction  technique  of  section  3.  The 
comparisons  will  be  made  directly  on  the  Green’s 
function,  since  this  is  the  relevant  physical  quantity 
rather  than  on  its  components,  the  Sommerfeld 
integrals.  From  the  asymptotic  evaluation  of  the 
Sommerfeld  integrals,  the  order  of  decay  of  each  of 
the  integrals  is  known  and  by  applying  this  informa¬ 
tion  to  the  Green’s  functions  of  (3)-(5),  the  surface 
wave  “radiation”  patterns  are  easily  obtained.  Fi¬ 
nally,  the  computation  times  of  the  numerical  and 
asymptotic  solutions  will  be  compared.  The  closed- 
form  asymptotic  solution  will  be  shown  to  be  ex¬ 
tremely  efficient  which  is  very  useful  for  the  numer¬ 
ical  computation  of  mutual  coupling  of  electrically 
large  slot  arrays  that  otherwise  might  have  been 
intractable. 

By  examining  the  asymptotic  solutions  given  by 
(3)-(5),  (31)-(36),  and  (43)^44),  one  can  determine 
that  the  dominant  field  contribution,  away  from  the 
source,  is  dictated  by  the  SWP  residue  terms  (that 
is,  surface  waves)  which  decay  as  The  SDP 
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contribution  decays  as  p~^  since  the  slope  transi¬ 
tion  function  approaches  unity  as  p  becomes  large. 
Because  of  the  dominance  of  the  surface  wave 
terms  (for  purposes  of  examining  the  far-field  “ra¬ 
diation”  pattern  of  a  magnetic  current  and  observa¬ 
tion  point  on  the  ground  plane),  it  is  only  necessary 
to  consider  the  residue  contributions.  This  implies 
that  the  Sommerfeld  integrals  of  type  1  are  approx¬ 
imately  the  negative  of  the  Sommerfeld  integrals  of 
type  2.  This  relationship  along  with  the  prescription 
of  constructing  the  Green’s  dyadic  in  (3)-(5)  yields 
the  copolarization  and  cross-polarization  surface 
wave  radiation  patterns  in  Figure  7.  The  patterns  in 
both  Figures  7a  and  7b  are  due  to  an  jc-oriented 
magnetic  current  source.  The  pattern  in  Figure  7a  is 


Figure  7.  Surface  wave  “radiation”  patterns  for  an 
AT-directed  magnetic  current  source  receiving  (a)  the  x-po- 
larized  field  (co-pol)  and  (b)  the  y-polarized  field  (cross- 
pol). 
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Fi^re  8.  Comparison  of  closed-form  uniform  asymp¬ 
totic  component  of  the  dyadic  Green’s  function  with 
“exact”  solution  computed  using  the  envelope  extraction 
technique.  Cut  taken  in  direction  of  surface  wave  null. 
The  parameters  e,  =  2.2,  /a,  =  1.0,  and  h  =  O.IAq. 


proportional  to  the  x  component  of  the  magnetic 
field,  and  the  pattern  in  Figure  lb  is  proportional  to 
the  y  component  of  the  magnetic  field.  Unlike  this 
sketch,  the  various  surface  waves  typically  have 
different  relative  magnitudes. 

The  first  example  of  a  closed-form  uniform  as¬ 
ymptotic  calculation  is  of  a  material  with  a  relative 
permittivity  of  2.2,  a  relative  permeability  of  1.0, 
and  a  dielectric  thickness  of  0.1  free  space  wave¬ 
length.  This  structure  has  a  single  TM  SWP  and  a 
single  TE  LWP  which  need  to  be  incorporated  into 
the  uniform  asymptotic  analysis  [Marin  et  al., 
1990].  Since  there  is  only  one  pole  in  the  proximity 
of  the  saddle  point  in  each  of  the  two  TM  integrands 
and  the  two  TE  integrands,  the  single  pole  formu¬ 
lation  of  section  4.1  was  used.  Only  one  surface 
wave  pole  exists  in  each  of  the  TM-type  integrals; 
hence  the  surface  wave  will  propagate  along  the 
dielectric  slab  broadside  to  the  magnetic  current. 
Figure  8  shows  the  real  and  imaginary  parts  of  the 
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Lateral  Separation  (Xq) 


Figure  9.  Comparison  of  closed-form  uniform  asymp¬ 
totic  component  of  the  dyadic  Green’s  function  with 
“exact”  solution  computed  using  the  envelope  extraction 
technique.  Cut  taken  in  direction  of  surface  wave  peak. 
The  parameters  =  2.2  ,  =  1.0,  and  h  =  0.1  Ag. 


component  of  the  dyadic  Green’s  function,  G, 
in  the  direction  of  the  surface  wave  null  (endfire). 
The  solid  lines  are  the  numerical  reference  solu¬ 
tions,  and  the  dotted  lines  are  the  closed-form 
asymptotic  solutions.  As  we  can  see,  the  agreement 
is  excellent  for  source  and  field  point  separations 
down  to  a  third  of  a  free  space  wavelength  or  so. 
Figure  9  also  shows  the  real  and  imaginary  parts  of 
the  Gxx  component  of  the  dyadic  Green’s  function, 
G,  for  the  same  structure,  but  this  is  for  the  broad¬ 
side  direction  along  which  the  lone  surface  wave 
propagates.  It  is  evident  that  the  agreement  is  again 
excellent,  and  the  closed-form  asymptotic  solution 
in  this  case  also  proves  accurate  down  to  approxi¬ 
mately  a  third  of  a  free  space  wavelength. 

The  next  case  under  consideration  has  a  relative 
permittivity  of  9.6,  a  relative  permeability  of  1.0, 
and  a  dielectric  thickness  of  0.05  free  space  wave¬ 
lengths.  As  in  the  previous  case  there  is  only  one 
pole  in  each  of  TM  and  TE  integrals  that  need  to  be 


considered  [Marin  et  al.,  1990].  Figures  10  and  11 
show  the  comparisons  between  the  numerical  and 
closed-form  uniform  asymptotic  solutions  of  the 
Gxx  component  of  the  dyadic  Green’s  function  for 
the  endfire  and  broadside  separations,  respectively. 
The  agreement  is  better  than  for  the  previous  set  of 
examples  in  that  there  appears  to  be  a  deviation 
between  the  two  solutions  for  source  and  field  point 
spacings  less  than  0.2  of  a  free  space  wavelength. 

The  third  and  final  case  shown  in  Figures  12  and 
13  is  one  where  there  are  two  poles  that  influence 
the  TE  SDP  evaluation  and  one  that  influences  the 
TM  SDP  evaluation.  The  communication  [Marin  et 
al.,  1990]  predicts  that  only  one  TE  pole  needs  to  be 
taken  into  account;  however,  we  found  that  two  TE 
poles  must  be  incorporated  into  the  analysis  to 
obtain  accurate  results.  The  relative  permittivity  is 
2.2,  the  relative  permeability  is  1.0,  and  the  material 
thickness  is  0.2  free  space  wavelengths.  By  com¬ 
parison  with  the  previous  two  cases  there  is  much 
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Figure  10.  Comparison  of  closed-form  uniform  asymp¬ 
totic  Gxx  component  of  the  dyadic  Green’s  function  with 
“exact”  solution  computed  using  the  envelope  extraction 
technique.  Cut  taken  in  direction  of  surface  wave  null. 
The  parameters  =  9.6,  =  1.0,  and  h  =  O.OSAg. 
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Stronger  interference  between  the  various  (space 
wave  and  pole  wave)  radiation  mechanisms.  The 
influence  of  these  mechanisms  is  more  evident  in 
this  case  because  two  TE  leaky  wave  poles  are 
close  to  the  saddle  point  in  the  spectral  domain. 
Once  again,  the  agreement  is  excellent  down  to  a 
lateral  source  and  field  point  separation  of  approx¬ 
imately  one  third  of  a  free  space  wavelength. 

The  major  advantage  of  the  closed-form  uniform 
asymptotic  evaluation  of  the  Green’s  function  is  the 
computational  time  savings  compared  to  numerical 
integration  techniques.  The  CPU  time  to  compute 
all  four  Sommerfeld  integrals  (i.e.,  all  components 
of  the  dyadic  Green’s  function)  is  tabulated  in  Table 
1  for  various  lateral  separations  of  the  source  and 
observation  points.  It  is  clear  that  the  closed-form 
uniform  asymptotic  representation  ofifers  signficiant 
computational  time  savings  compared  to  numerical 
integration.  As  the  lateral  separation  increases,  the 
computation  time  of  the  numerical  integration 


Figure  11.  Comparison  of  closed-form  uniform  asymp¬ 
totic  Gjj  component  of  the  dyadic  Green’s  function  with 
“exact”  solution  computed  using  the  envelope  extraction 
technique.  Cut  taken  in  direction  of  surface  wave  peak. 
The  parameters  =  9.6,  =  1.0,  and  /i  =  O.OSAq. 
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Figure  12.  Comparison  of  closed-form  uniform  asymp¬ 
totic  component  of  the  dyadic  Green’s  function  with 
“exact”  solution  computed  using  the  envelope  extraction 
technique.  Cut  taken  in  direction  of  surface  wave  null. 
The  parameters  e,  =  2.2,  =  1.0,  and  h  =  0.2Ao. 


grows  without  bound,  whereas  the  computation 
time  of  the  asymptotic  solution  remains  constant. 
This  along  with  the  high  degree  of  accuracy  makes 
the  closed-form  uniform  asymptotic  solution  ideal 
for  the  analysis  of  not  only  large  slotted  material 
coated  planar  structures  but  also  for  relatively  small 
structures. 

In  summary,  two  methods  of  computing  the  dy¬ 
adic  aperture  Green’s  function  G  of  (2)  were  pre¬ 
sented.  One  was  a  numerical  evaluation  of  the  exact 
integral  representation  which  is  valid  for  all  values 
of  lateral  separation  of  the  source  and  observation 
points.  This  solution  incorporates  the  envelope  and 
singularity  extraction  techniques  to  enhance  the 
computation  efficiency.  The  second  solution  and 
the  major  contribution  of  this  paper  is  a  closed-form 
uniform  asymptotic  representation  that  is  valid  for 
large  lateral  source  and  field  point  separations; 
however,  it  was  shown  to  remain  valid  for  lateral 
separations  down  to  a  half  of  a  free  space  wave- 
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Figure  13.  Comparison  of  closed-form  uniform  asymp¬ 
totic  Gjiji  component  of  the  dyadic  Green’s  function  with 
“exact”  solution  computed  using  the  envelope  extraction 
technique.  Cut  taken  in  direction  of  surface  wave  peak. 
The  parameters  =  2.2,  =  1-0,  and  h  =  O.ZAq. 

length  and  even  two  tenths  of  a  free  space  wave¬ 
length  in  some  cases.  This  asymptotic  solution  is  at 
least  2  to  3  orders  of  magnitude  computationally 
more  efficient  than  the  numerical  representation  for 
distances  characteristic  of  large  arrays.  The  combi¬ 
nation  of  these  two  solutions  is  ideal  for  the  efficient 
calculation  of  the  aperture  dyadic  Green’s  function. 


Table  1.  CPU  Time  (Milliseconds  on  a  VAX 
8550)  to  Calculate  All  Elements  of  the  Aperture 
Dyadic  Green’s  Function  for  the  Case  of  =  9.6, 
fir  =  1.0  and  h  =  0.05Ao 


Lateral 

Distance, 

p/Ao 

CPU  Time 
Numerical 
Integration 

CPU  Time 
Closed-Form 
Asymptotic 

0.5 

141 

0.84 

1.0 

228 

0.84 

2.0 

335 

0.84 

10.0 

991 

0.84 

G,  for  the  entire  range  of  lateral  separation.  The 
results  developed  here  can  be  used  for  a  variety  of 
slot/feed  geometries  containing  material  layers  be¬ 
low  (substrate)  and  a  material  layer  above  (super- 
strate)  the  slots  in  a  PEC  ground  plane. 
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EM  PLANE-WAVE  DIFFRACTION  BY  A 
MATERIAL-COATED  PERFECTLY 
CONDUCTING  HALF-PLANE<- 
OBLIQUE  INCIDENCE 
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Department  of  Electrical  Engineering 
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Columbus,  Ohio  43212 
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Coaled  half-plane,  diffraction,  Wiener-Hopf,  generalized  impedance 
boundary  conditions 

ABSTRACT 

A  high-frequency  solution  is  developed  for  the  electromagnetic  dif¬ 
fraction  of  an  obliquely  incident  plane  wave  by  a  material-coated 
perfectly  conducting  half-plane.  It  is  assumed  that  the  coating  is 
electrically  thin  so  that  it  can  be  modeled  by  a  generalized  imped¬ 
ance  boundary  condition  of  0(t),  where  t  is  the  coating's  thickness. 
©  1994  John  Wiley  i  Sons,  Inc. 

I.  INTRODUCTION 

The  canonical  problem  to  be  studied  is  the  high-frequency 
electromagnetic  (EM)  diffraction  by  a  material-coated  per- 
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(d)  (e) 


Figure  1  Material-coated  PEC  half-plane  geometry  and  equivalent  configurations 


fectly  conducting  (PEC)  half-plane,  with  the  same  coatings 
on  the  top  face  {jr  >  0,  y  =  0*}  and  the  bottom  face  {x  >  0, 
>’  =  0  ).  The  incident  field  is  a  plane  wave  of  arbitrary  po¬ 
larization  obliquely  incident  to  the  z  axis  with  an  angle  %' 
(0  <  e'  <  rr).  The  geometry  pertaining  to  this  problem  is 
shown  in  Figure  1(a),  with  <t>  and  <t>'  being  the  scattered  and 
incident  angles,  respectively.  It  is  important  to  keep  in  mind 
that  the  half-plane  being  considered  here  is  opaque,  whereas 
the  one  considered  in  Section  V  of  [1]  is  transparent. 

For  the  coated  PEC  half-plane  depicted  in  Figure  1(a)  with 
the  same  material  coating  on  both  faces  {y  =  0=},  the  ge¬ 
ometry  is  symmetrical  about  the  y  =  0  plane.  Hence,  the 
coated  half-plane  may  be  decomposed  into  two  equivalent 
configurations  by  first  expressing  the  incident  field  (£',  H'^) 
as  a  superposition  of  the  even  and  odd  symmetrical  compo¬ 
nents  (Figures  1(b)  and  1(c)].  These  configurations  are  in  turn 
equivalent  (for  y  >  0)  to  the  ones  shown  in  Figures  1(d)  and 
1(e);  a  PEC  ground  plane  with  semi-infinite  coating  and  a 
coated  PEC  half-plane  joined  with  a  PMC  (perfectly  magnetic 
conducting)  half-plane.  Note  that  the  latter  geometry  is  not 
considered  in  (!].  The  solution  to  the  geometry  in  Figure  1(a) 
may  then  be  obtained  via  a  superposition  of  the  results  for 
Figures  1(d)  and  1(e). 

The  coating's  thickness  t  is  assumed  to  be  electrically  small, 
so  that  the  employment  of  a  generalized  impedance  boundary 
condition  (GIBC)  of  0(i)  provides  an  accurate  model  for  the 
coated  half-plane.  This  GIBC  of  0(r)  may  be  expressed  as 


follows  (n  =  ±^1,  y  S  0)  [2]: 

=  -ikrjol  n  X  //(y  =  0=),  x  >  0,  (1) 

where  170  and  k  are,  respectively,  the  free-space  impedance 
and  wave  number.  Note  that  the  GIBC  of  0(/)  given  above 
can  be  rewritten  in  terms  of  the  z-component  fields  £.  and 
//j  only,  because  the  planar  geometry  in  Figure  1(a)  is  two 
dimensional,  though  the  incidence  is  oblique.  Also,  it  has  been 
shown  by  Ly,  Rojas,  and  Pathak  [1]  that  a  complete  solution 
for  a  boundary-value  problem  involving  GIBC  of  0(/)  can  be 
obtained  by  imposing  a  junction  (or  contact)  condition  which 
may  be  found  via  a  quasistatic  analysis.  This  junction  con¬ 
dition  is  applied  in  addition  to  the  well-known  edge  condition. 
The  solution  thus  obtained  automatically  satisfies  reciprocity 
even  for  the  case  of  oblique  incidence. 

II.  PEC  GROUND  PLANE  VITITH  A  SEMI-INFINITE 
MATERIAL  COATING 

The  partially  coated  ground  plane  shown  in  Figure  1(d)  is  a 
special  case  of  the  more  general  two-part  coated  ground  plane 
discussed  by  Ly  et  al.  [1],  The  configuration  in  Figure  1(d)  is 
obtained  from  the  one  depicted  in  Figure  2  in  [1]  as  the  thick- 
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ness  of  the  material  coating  {(2)  on  the  half-plane  {y  =  0, 
X  <  0}  shrinks  to  zero.  Thus,  by  taking  appropriate  limits  of 
the  corresponding  parameters  in  the  solution  for  the  two-part 
coated  PEC  ground  plane  analyzed  in  [1],  one  may  obtain 
the  spectral  function  .'4(h',  <I>';  6')  (i.e.,  the  spectral  scattered 
field  in  the  w  domain  that  would  yield  a  spatial  field  via  an 
inverse  Fourier  transform)  for  the  partially  coated  ground 
plane  as  follows; 


■JXw,  4>'-,  0')  =  C(k-)T(w)  sin  w 

^  [cos  w  +  cos  4)' 

+  5(w,  <t>'-  e')j  sin  4>'T(<#.’)C(<^')Fo„  (2) 


where  ^02  -  ^02  and  H02  are,  respectively,  the  mag¬ 

nitudes  of  the  incident  fields  £j  and  H'^  at  the  origin.  The 
matrices  in  (2)  are  defined  as  follows: 


maining  terms  present  in  the  above  expressions  are  given  by 


sin  C  = 


-  l)t’ 
a  =  sin  -I-  sin  Vj,  = 


sin  Po  - 

b  2 


-  b,  (5a) 


Ifr 


Kie,  -  1)/’ 


r  Mrfr  -  1 

b  =  sm  Pa  sin  Pi,  =  — :  ,  ,,, - 77. 

Sin^  e  (€r  -  1) 


(5b) 


III.  PLANAR  JUNCTION  OF  COATED  PEC  HALF-PLANE 
WITH  PMC  HALF-PLANE 

The  geometry  pertaining  to  this  scattering  problem  is  illus¬ 
trated  in  Figure  1(e),  where  one  half-plane  is  a  PMC  and  the 
other  is  a  coated  PEC.  Although  this  geometry  is  not  consid¬ 
ered  in  [1],  it  can  still  be  obtained  from  the  two-part  material- 
coated  PEC  ground-plane  geometry  of  Figure  2  in  [1]  by 
means  of  a  key  step.  In  other  words,  the  coated  PEC  half¬ 
plane  {x  <  0,  y  =  0}  in  Figure  2  in  [1]  becomes  a  PMC  half- 


,  1  -  cos  w  sin  w  cos  6  . ,  ,  ?  ,  7  •  ?  /,- 

Vi')  =  ,  Afw)  =  cos^  6  +  cos^  w  sin'^  6  ; 

A(vi’)  sin  w  cos  C  cos  w  '  ' 


C( 


^(w)  = 


0 


0  g-(H',  »'a)g,(vt',  Pk) /[Vicos  J 


A(w,  0')  =  °  1,  S(W,  4>';  0')  =  b"']  ; 

U  A22}  "22. 


in  which 


All  =  0'  cos  vv  cos  <f>'  -  cos^  6', 

A22  =  (cos  w  cos  4>'  +  cos^  wo)(cos^  €'  -  sin^  6'  cos  w  cos  <^') 


,  r  ■  7  ,,  cos  Wo  ,  /cos  w  cos  <f)' 

+  (cos  w  +  cos  <p  )  sin^  0  cos  w  cos  <f>  — — r-  +  cos^  0  I  — : - 

L  Kwo) 


cos  w  -  cos  <t>‘ 


Bn  = 

Bi2  = 
B21  = 
B22  = 

A(**'o)  = 
i‘(wi.)  = 


1  cos  g'[ie(w;0  -t-  i^(w,')]  -t-  /i(wo)  cos  Wq  sin  g'[iXw;-)  -  2e(w,~)]  . 

/i(wo)  cos  Wq  sin  fl'[i’(w/')  +  i^(w,‘')]  +  1  cos  fl'{y(w;.)  -  if(w;.)] 
sin  0'  A(wo)(l/sin^  6'  +  fl/sin  6’  -I-  b)"' 
b(wo)  cos  Wo  sin  ^'[^(w;.)  -t-  ie(w,')]  -f  i  cos  ®'[iXw;-)  -  Jf(w;.)] 

_ sin^  fl'(l/sin  0"  sin  _ 

b(wo)  cos  Wo  sin  e’[i?(H'^)  +  ^( *♦'«  )]  +  <  cos  fl'(Jf(w;-)  - 

_ if(w;.)  -  i?(w.-) _ 

h(wo)  cos  Wo  sin  b'[i’(**'^)  +  ^( **'«•)]  +  i  cos  ®'(y(w;')  -  **'»  )] 

•  (z  cos  6'  A(wo)/(A(wo)  cos  Wo)][/i(wo)  cos  Wo  -  cos  <b'][b(wo)  cos  wo  -  cos  w]  (4f) 
fl(l  -  cos  Wo)g!(Wo,  »'„)g^,(Wo,  Pb)  +  1 


\h{wo)  cos  Wo 

1  cos  0'  sin  0', 

2  cos  6'[/l(Wo)  cos  Wo  -  cos  <b’], 
2  cos  P'[b(wo)  cos  Wo  -  cos  w]. 


■) 


(3a) 

(3b) 

(3c) 

(4a) 

(4b) 

(4c) 

(4d) 

(4e) 


fl(l  -  cos  Wo)gi(Wo,  l<,)gU*^0,  h)  -  1’ 
g-(>^^  Va)gA^i',  Vb)g-(*vi',  0 
V2  cos(w^,/2) 


(4g) 

(4h) 


Note  that  the  expressions  for  the  functions  g-  were  summa¬ 
rized  in  [3]  and  are  not  repeated  here.  Also,  K  =  k  sin  0", 
cos^  Wo  =  1  +  b,  wi'  =  ir/2  i  i  ln(cot(e72)),  and  the  re¬ 


plane  without  any  coating  by  applying  the  limits  -♦  0,  and 
(Mfj  ~  l)f2  ^  “  (which  is  the  key  step  in  this  analysis)  to  the 
solution  for  the  two-part  material-coated  PEC  ground  plane 
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IV.  FINAL  RESULT 

A  superposition  of  the  spectral  functions  If  and  If  given  in  the 
previous  sections  yields  the  spectral  field  for  the  material- 
coated  PEC  half-plane  with  the  same  coatings  on  both  faces. 
Using  the  steepest  descent  method  to  evaluate  the  inverse 
Fourier  transform  of  the  spectral  field,  one  may  express  the 
total  field  for  the  coated  PEC  half-plane  as  a  sum  of  the 
inddent-,  reflected-,  diffracted-,  and  surface-wave  fields.  The 
diffracted  field  at  the  point  (p,  0,  w  -  ff')  may  be  written  as 


^5,(101,0';  O')!  +  2:?. 

0  J  i-l 


.  1  -  F(-2Kpy^) 

2r 

0  <  |0|,  0'  <  IT, 


]■ 


(8) 


where 


-  [o  10!/ 0]’  “"  (  '2 


i  =  1.  2,  3, 


(9) 


(b)  TEji  poiftriMtiofi:  £oj  s  0^  «  1. 

Figure  2  Total  field  for  a  material-coated  PEC  half-plane  excited 
by  an  obliquely  incident  plane  wave,  (a)  TM,  polarization:  £o,  =  1, 
ZuHit-  =  0.  (b)  TE,  polarization:  Eo,  =  0,  •=  1 

alluded  to  above.  The  above  limits  yield  the  following  spectral 
field  for  the  present  problem; 

7f.(H',  0';  ff)  =  C(h’)^'(**')  sin  **'  [  ^  ^ 

^  \  \  w  +  cos  4> 

+  ff'fH',  0';  0')]  •  sin  0'T'(0')ZC(0’)Fo„ 
where  the  new  matrices  are  defined  by 


(6) 


and  Wi  =  IT  -  0',  Wj  =  -vj,,  *V3  =  The  functions 
r,  i  =  1,  2,  3,  are  the  residues  of  the  spectral  function 
L(if;(H',  0';  O')  +  |0l/0if,(»v,  0';  »')]  at  the  poles  w  =  Wj. 
The  residue  Fj  also  produces  a  surface-wave  field  given  by 

_  101),  (10) 

where  0,  =  -Re  -I-  arccos[l/cosh(Im  v^)].  The  geomet¬ 
rical  optics  fields  are  not  included  here  due  to  space  limita¬ 
tions.  A  numerical  example  based  on  the  high-frequency  so¬ 
lution  obtained  above  is  shown  in  Figure  2,  where  the  plane 
wave  is  obliquely  incident  {O'  =  45°)  to  the  edge  of  the  coated 
half-plane.  Both  fields  £,  and  H,  are  excited  by  the  incident 
field  of  either  TM,  or  TE,  polarization,  though  the  cross- 
polarized  field  appears  small  for  most  scattering  aspects.  The 
effect  of  the  surface  wave  can  clearly  be  seen  in  Figure  2 
around  0  =  360®,  where  the  total  field  starts  to  increase  in 


AXw,  0';  e')  -  5'(*v,  0'; 


All  =  sin^  g'  cos  w  cos  0'  -  cos^  O',  =  a  sin^  0’  cos  w  cos  0'  -  cos^  O', 

r{wt)  -  r(w;) 


B'li  =  -i  cos  ff  sin  ff 


B'n  =  2a  cos  9'  sin  0' 


B21  =  -  2  cos  9"  sin^  9' 


-h  s'iw^Y 

,  (1/sin^  9*  +  g/sin  O'  -f  b)~* 


£'(**'?)  +  f'(w;.) 

(1/sin  9*  -I-  sin  ()■’ 

f'(H-;.)  +  f'(w;.)  ’ 


=  ia  cos  9'  sin  9' 


f'(w;)  -I-  i-XH-;')’ 


r{wp)  = 


y/2  sin(»’J-/2) 


(7a) 

(7b) 

(7c) 

(7d) 

(.7e) 

(7f) 

(7g) 

(7h) 
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magnitude  due  to  the  surface  wave  field,  especially  for  the 
case  of  larger  e,  and  fi,.  Note  that  as  the  oblique  incidence 
becomes  normal  incidence  (i.e. ,  ff  -*  90°),  the  result  obtained 
here  reduces  to  the  Maliuzhinets-based  solution  developed 
by  Ly  and  Rojas  [4]  when  the  expressions  in  [4]  are  simplified 
to  the  special  case  of  equal  coatings  on  both  faces  of  the  half¬ 
plane.  The  accuracy  of  the  two-part  material  slab  solution 
from  which  the  present  solution  is  obtained  has  been  verified 
in  (1)  using  a  dependent  moment-method  solution  as  a  ref¬ 
erence. 
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Integra]  Equation  Analysis  of  a  Sheet  Impedance 
Coated  Window  Slot  Antenna 

R.  Torres  and  E.  H.  Newman 


A^sinicl— This  paper  will  present  an  integral  equation  and  method  of 
moments  analysis  of  a  window  slot  antenna.  The  antenna  is  modeled  by  a 
sheet  admittance  coated  rectangular  aperture  in  an  infinite  ground  plane. 
It  is  shown  that  the  sheet  admittance  coated  aperture  is  complementary 
to  a  sheet  impedance  plate,  which  permits  the  window  slot  antenna 
to  be  analyzed  with  existing  computer  programs.  Numerical  results 
are  compared  with  measurements  for  input  impedance  and  radiation 
efficiency. 


I.  Introduction 

This  paper  will  present  an  integral  equation  and  method  of  mo¬ 
ments  (MM)  [1]  solution  for  a  recently  developed  rectangular  slot 
antenna,  which  can  be  mounted  in  the  outer  perimeter  of  a  window 

[2] ,  f3].  A  conventional  slot  antenna  consists  of  a  slot  or  gap  cut 
in  a  perfectly  conducting  ground  plane  [4],  [5].  If  the  cross-hatched 
regions  of  Fig.  1  are  considered  to  be  perfectly  conducting  surfaces, 
then  Fig.  1  represents  a  conventional  rectangular  slot  antenna  placed 
in  the  outer  perimeter  of  a  rectangular  window.  The  problem  with 
this  design  is  that  it  is  not  possible  to  sec  through  the  window  since 
it  has  been  covered  with  a  perfect  conductor.  Recently,  Torres  and 
Walton  have  solved  this  problem  by  replacing  the  irmer  perfectly 
conducting  patch  by  a  resistive  surface  of  typically  4-8  d/D  [2], 

[3] .  These  resistive  surfaces  can  be  made  almost  transparent  at 
optical  ftequencies,  and  yet  have  low  enough  resistance  to  serve  as 
a  reasonable  ground  plane  for  the  slot  antenna.  For  example,  Torres 
and  Walton  designed  an  automobile  windshield  slot  antenna  which  is 
typically  5  dB  below  the  standard  whip  in  the  FM  band,  but  about 
3  dB  above  the  whip  in  the  AM  band. 

The  next  section  develops  an  integral  equation  for  a  sheet  admit¬ 
tance  coated  aperture  in  an  infinite  ground  plane,  which  serves  as 
an  idealized  model  for  the  window  slot  antenna  since  it  neglects 
antenna/automobile  interactions.  Next,  it  is  shown  that  the  sheet 
admittance  coated  aperture  is  complementary  to  a  sheet  impedance 
plate,  which  allows  the  window  antenna  to  be  analyzed  with  existing 
computer  codes  [6].  Finally,  numerical  results  are  compared  with 
measurements  for  the  input  impedance  and  radiation  efficiency  of  the 
window  slot  antenna. 


n.  TReory 


A.  Iruroduction 

This  section  will  develop  an  integral  equation  for  a  sheet  impedance 
coated  slot  antenna.  As  illustrated  in  Fig.  1,  the  antenna  consists  of  a 
rectangular  aperture  of  width  W  and  height  I  in  an  infinite  perfectly 
conducting  ground  plane.  The  outer  edge  of  the  rectangular  ^lerture 
contains  a  gap  of  width  g  and  of  sheet  impedance  Z,g,  while  the 
remaining  center  section  (window)  of  the  aperture  contains  the  sheet 
impedance  Z,c.  The  term  sheet  impedance  implies  a  generalization 
of  a  resistive  surface,  where  the  resistivity  may  be  complex. 

Manuscript  received  September  8,  1992;  revised  June  9,  1993. 
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Window  Aritenno 


Fig.  1.  Geometry  for  a  sheet  impedance  coated  window  antenna. 


A  sheet  impedance  is  a  zero  thickness  model  for  an  electrically  thin 
dielectric  slab  [7]-[]4].  In  particular,  an  electrically  thin  dielectric 
slab  of  thickness  T,  relative  permittivity,  fr,  conductivity  a,  and 
complex  permittivity  e  =  Crfo  —  j{a/u))  has  sheet  impedance 


Z.= 


(n/o) 


(1) 


where  «o  is  the  permittivity  of  free  space  and  w  is  the  radian 
frequency.  Note  that  for  a  lossless  dielectric,  e  is  pure  real,  and 
thus  Z,  is  pure  imaginary.  By  contrast,  for  a  very  lossy  dielectric, 
<r/w  >  <rfo,  and  Z,  =  1/(tT  is  pure  real.  A  pure  real  Z,  is 
commonly  referred  to  as  a  resistive  surface  or  a  space  cloth.  A  surface 
of  zero  sheet  impedance  is  a  perfectly  conducting  surface,  while  an 
infinite  sheet  impedance  corresponds  to  the  absence  of  the  dielectric 
slab.  If  Z,g  =  oo  and  Z,c  =  0,  then  the  antenna  of  Fig.  1  would 
be  a  conventional  rectangular  slot  antenna.  The  volume  equivalence 
[15]  can  be  used  to  show  that  the  sheet  inqiedance  surface  can  be 
replaced  by  free  space  and  the  equivalent  electric  surface  current 


J  =  E/Z.  (2) 

where  E  is  the  electric  field  at  the  surface  of  the  sheet  impedance.  It 
is  implicit  that  (2)  ^lies  only  to  the  components  of  E  and  J  which 
are  tangential  to  the  sheet  impedance  surface. 

For  multilayer  dielectrics,  the  individual  sheet  impedance  of  the 
layers  ^>pear  as  impedances  in  parallel.  For  example,  the  sheet 
admittance  of  a  two-layer  electrically  thin  dielectric  slab  is  simply 
[8] 

y-  =  =  Y*i+Y*2  =  Mfi-fo)Ti+M€3~(o)T2  (U/D)  (3) 

where  tj  and  Ti  or  <2  and  T2  are  the  permittivity  and  thickness  of 
layer  I  or  2,  respectively.  Equation  (3)  can  be  used  to  determine  the 
sheet  inqiedance  of  a  thin  plastic  window  coated  with  a  resistive  film. 


B.  Sheet  Admittance  Apei^re  Integral  Equation 

Fig.  2(a)  shows  the  side  view  of  a  planar  ^lerture  A  in  an 
infinite  and  perfectly  conducting  ground  plane.  The  regions  above 
and  below  the  ground  plane  will  be  referred  to  as  Regions  a  ind  b, 
respectively.  An  infinitesimal  distance  6  above  or  below  the  aperture 
are  the  nonuniform  sheet  admittances  y,«  and  Y,k,  respectively.  The 
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free  space  equivalent  problenis  of  Fig.  2(c)  and  (d). 

E'  +  B'”'>  +  H{2M)  +  {ff[Jo(at  y  =  ^)] 

+  flr[-Ja(aty  =  -«)]} 

=  F(-2Af)  +  {F[-/6(at  y  =  «)] 

+  ff  [J6(at  y  =  -^)]}  (tang,  in  A)  (7) 

where  is  the  free  space  magnetic  field  of  the  image  impressed 
currents.  At  almost  every  point  in  space,  the  electromagnetic  fields 
of  Jo(aty  =  6)  are  exactly  cancelled  by  the  image  currents 
— Ja(at  y  =  —6).  The  one  exception  is  for  field  points  in  the  aperture 
A  where  boundary  conditions  can  be  used  to  show  that  the  free  space 
tangential  magnetic  fields  of  Jo(at  y  =  6)  and  -Ja(at  y  —  —6)  are 
identical  and  equal  to  j  x  Jo/2.  Then,  using  (2)  and  (4), 

fl-(Jo(at  y  =  6)]  +  F[-Jo(at  y  =  -«)] 

=  y  X  Jo  =  »  X  y.o£  =  -Y.aM.  (8) 


Similarly, 


(e)  (d) 

Fig.  2.  (a)  An  aperture  coated  with  a  sheet  admittance,  (b)  and  (c)  Equivalent 
problems  in  Region  a.  (d)  An  equivalent  problem  in  Region  6. 

excitation  is  by  the  impressed  currents  (J*,  M')  which  are  assumed 
to  be  in  Region  a.  (J',  M')  produce  the  unknown  total  fields  {E,  H) 
in  either  Regions  a  or  6  when  radiating  in  the  presence  of  the  aperture 
and  the  sheet  admittances,  and  the  assumed  known  incident  fields 
{E'.  H')  in  free  space.  In  this  pjqier,  all  fields  and  currents  are  time 
harmonic,  with  the  e^“‘  time  dependence  suppressed. 

Following  the  usual  procedure  [16],  the  first  step  in  deriving  the 
aperture  integral  equation  is  to  construct  equivalent  problems  in 
Regions  a  and  b.  An  equivalent  problem  for  Region  a  is  shown  in  Fig. 
2(b)  in  which  the  aperture  is  shorted  and  the  magnetic  surface  current 

.  M  =  E  X  y  in  A  (4) 

is  placed  just  above  the  aperture.  Also,  (2)  is  used  to  replace  the 
sheet  admittance  l.a  by  free  space  and  the  equivalent  electric  surface 
current 

Ja  =  Y,eE  (comp.  tang,  to  A).  (5) 


JBr(-J6(at  y  =  «)]  +  F[J6(at  y  =  -6)] 

=  — y  X  Ji  =  — y  X  Y.bE  =  Y,hM,  (9) 

and  where  it  is  emphasized  that  (8)  and  (9)  apply  to  the  tangential 
components  in  A.  Inserting  (8)  and  (9)  into  (7),  and  using  the  fact 
that  in  the  plane  of  the  aperture  the  tangential  components  of  W  and 
E'”"^  arc  equal  yields 

2^(2Af)  -  (y„  +  y,6)M  =  -2JT’  (tang,  in  A).  (10) 

As  expected,  M  does  not  depend  upon  how  the  sheet  admittance 
is  distributed  above  or  below  the  aperture,  but  instead  is  dependent 
upon  the  sum  or  parallel  combination  of  the  sheet  admittances  yjo 
and  Y,k.  Defining  the  net  sheet  admittance  as 

y.  =  y.a+v,t,  (11) 

(10)  becomes 

B{2M)  -  ^(2M)  =  -B'  (tang,  in  A).  (12) 

4 

Equation  (12)  is  the  desired  integral  equation  for  the  equivalent 
magnetic  surface  current  representing  the  sheet  admittance  coated 
aperture.  Note  that  it  has  been  written  in  terms  of  2M  rather  than  M 
since  the  scattered  fields  are  the  free  space  fields  of  +2Af  or  —2M 
in  Regions  a  or  6,  respectively. 


In  the  equivalent  problem  of  Fig.  2(b),  the  superposition  of  the  fields 
of  (J’,  Af‘),  3f,  and  Ja,  radiating  in  the  presence  of  the  infinite 
perfectly  conducting  ground  plane,  produce  the  proper  total  fields 
{E.  B)  above  the  ground  plane  and  a  null  field  below  the  ground 
plane.  In  the  free  space  equivalent  problem  of  Fig.  2(c),  image  theory 
has  been  used  to  remove  the  ground  plane  and  replace  it  with  the 
images  of  the  above  currents.  In  Region  a  of  the  free  space  equivalent 
problem  of  Fig.  2(c),  the  total  fields  {E,B)  are  produced  by  the 
superposition  of  the  free  space  fields  of  the  currents  (J',  M'),  and 
Af,  and  Ja  plus  their  images.  Following  the  same  procedure.  Fig. 
2(d)  shows  the  free  space  equivalent  problem  in  Region  b.  In  Region 
6,  the  total  fields  are  produced  by  the  sup)erpx)sition  of  the  free  space 
fields  of  the  magnetic  surface  current  —2M  located  in  A,  the  electric 
surface  current 

J(,  =  Y,i,E  (comp.  tang,  to  A)  (6) 

located  a  distance  S  below  the  aperture,  and  the  current  —  Js  located 
a  distance  6  above  the  apierture. 

The  aperture  integral  equation  is  an  expiression  of  continuity  of  the 
total  tangential  magnetic  field  across  the  apierture.  Referring  to  the 


C.  Complementary  Sheet  Impedance  Plate 
This  section  will  show  that  the  integral  equation  for  the  sheet 
admittance  coated  aperture,  (12),  is  complementary  to  that  of  a  sheet 
impedance  plate.  This  comparison  is  useful  since  it  is  a  relatively 
simple  matter  to  modify  one  of  the  many  existing  MM  computers 
for  perfectly  conducting  plates  so  that  they  can  treat  sheet  impedance 
plates  [16]. 

Consider  electric  and  magnetic  currents  related  by  J'  =  M'/r). 
Directly  from  Maxwell’s  equations  it  can  be  shown  that  the  elec¬ 
tromagnetic  fields  of  J'  and  M',  in  a  homogeneous  medium  with 
characteristic  impedance  t),  are  simply  related  by 

E(J')  -  7iS(M')  and  B{J')  =  -E{M')lr,.  (13) 

Eiefining  the  complementary  electric  surface  current  density  as  Jc  = 
-2M/ri,  using  (13),  and  multiplying  by  »?,  (12)  becomes 

-E{Jc)  +  ^Jc  =  -vB'  (tang,  in  .4).  (14) 

Now,  consider  a  sheet  impedance  plate,  with  sheet  impedance  Z.p  on 
the  plate  surface  S,  and  excited  by  the  impressed  currents  ( Jj,,  AfJ,), 
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with  free  space  fields  (iEj,,  E),).  Jp  is  given  by  the  sheet  impedance 
integral  equation  [7H14] 

~E(,Jp)  +  Z.pJp  =  E'p  (tang,  on  S).  (15) 

Assuming  that  the  impressed  cunents  for  the  sheet  impedance  plate 
are  the  complement  of  those  for  the  sheet  admittance  iperture,  i.e., 
(Jp,  M'p)  —  (-M'/ti,  then  the  incident  fields  will  be  related 
by 

(E'p,  E'p)  =  (-f,E\  E'h).  (16) 

In  this  case,  (15)  becomes 

-E(Jp)  +  Z,pJp  =  -riE'  (tang,  on  S),  (17) 

Comparing  (14)  and  (17)  shows  that  if  the  tperture  area  A  is 
identical  to  the  plate  surface  5  and 

2 

Z,p($bea  impedance  plate)  =  —y, (sheet  admittance  aperture), 

(18) 

then 

Jp  =  Jc  =  -2M/7}.  (19) 

That  is,  the  equivalent  magnetic  surface  current  2M,  representing 
the  sheet  admittance  aperture,  is  the  complement  of  the  equivalent 
electric  surface  current  Jp,  representing  the  complementary  sheet 
impedance  plate. 

in.  Numerical  Results 

This  section  will  present  numerical  results  for  the  input  impedance 
and  radiation  efficiency  of  the  window  slot  antenna.  All  numerical 
results  were  obtained  by  solving  the  complementary  equation  (15) 
with  the  authors’  Electromagnetic  Surface  Patch  Code:  ESP  IV  [6]. 
Referring  to  Fig.  1,  all  data  are  for  a  window  slot  antenna  of 
dimensions  =  £  =  31  mm  and  with  an  air  gap  (Z,p  =  oo) 
of  width  5  =  1  mm.  For  the  measured  results,  the  infinite  perfectly 
conducting  ground  plane  was  simulated  wiA  a  90  cm*  aluminum 
plate. 

Fig.  3  shows  a  comparison  of  the  computed  (solid  line)  and 
measured  (dashed  line)  input  impedance  of  the  window  slot  antenna 
for  inner  patches  of  resistivity  Z,c  =  0, 25  n/D,  and  200  Q/D  from 
2500  to  3200  MHz.  Over  this  frequency  range,  the  total  gap  length 
varies  from  lA  to  1.28A. 

Fig.  4  shows  a  comparison  of  computed  and  measured  (using 
the  Wheeler  cap  method  [17])  radiation  efficiency  as  a  function 
of  the  resistivity  of  the  inner  patch  at  /  :=  2.5  GHz.  Measured 
values  were  obtained  only  for  Z,c  =  25  fl/D  and  200  n/O,  and 
are  in  reasonable  agreement  with  computations.  As  expected,  as  the 
resistivity  of  the  inner  patch  increases  from  the  zero  resistance  of 
a  perfectly  conducting  surface,  the  radiation  efficiency  decreases. 
However,  for  very  large  resistivity,  the  resistance  patch  is  essentially 
approaching  free  space,  and  the  efficiency  begins  to  increase. 

rv.  Summary 

This  paper  has  presented  an  integral  equation  and  method  of  mo¬ 
ments  solution  for  a  window  slot  antenna.  The  window  slot  antenna 
is  modeled  by  a  rectangular  aperture  with  a  sheet  admittance  coating. 

It  is  shown  that  the  sheet  admittance  ipeiture  is  complementary  to  a 
sheet  impedance  plate,  which  permits  the  window  slot  antenna  to  be 
analyzed  with  existing  computer  codes  [6].  Numerical  results  for  the 
input  impedance  and  radiation  efficiency  are  found  to  be  in  reasonable 
agreement  with  measurements. 


Frequency  (MHz) 


Fig.  3.  Comparison  of  measured  and  computed  input  impedance  versus 
frequency  for  a  window  slot  antenna  with  resistivity  Z,c  =  0.  250/0, 
and  200  0/0  for  the  inner  patch. 


Fig.  4.  Comparison  of  measured  and  computed  radiation  efficiency  versus 
the  resistivity  of  the  inner  patch  at  /  =  2.5  GHz. 
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Integra]  Equation  Analysis  of  Artificia]  Media 

J.  L  Blanchard,  E.  H.  Newman,  and  M.  E.  Peters 


Abanet — This  paper  presents  an  integral  equation  and  method  of 
moments  (MM)  solution  to  determine  the  effective  permittivity  and 
permeability  of  an  artificial  medium.  The  artificial  medium  is  modeled 
as  a  triple  infinite  periodic  array  of  identical  scattering  elements.  A 
plane  wave  of  unknown  phase  constant  is  assumed  to  propagate  in  the 
periodic  medium  in  a  given  direction,  and  the  periodic  aaoment  method 
(PMM)  is  used  to  set  up  a  matrix  equation  for  the  currents  on  the  eroter 
clement  of  the  periodic  array.  By  setting  the  determinant  of  the  PMM 
impedance  matrix  to  sero,  one  can  determine  the  phase  constant  of  the 
|dane  wave,  and  then  the  effective  permittivity  and  permeability  ti  the 
artificial  m^um. 


I.  iNTllODUCnON 

An  artificial  medium  consists  of  a  large  number  of  scattering 
elements  placed  in  some  homogeneous  background  or  host  medium 
tl].  Typically  there  will  be  a  large  number  of  identical  electrically 
small  scatteren  per  cubic  wavelength.  When  an  electromagnetic  wave 
propagates  through  the  artificial  medium,  it  will  induce  currents 
on  or  in  the  scatterers.  Thus,  the  scatterers  can  be  viewed  as 
macroscopic  electric  and/or  magnetic  dipole  moments,  in  analogy  to 
die  microscopic  molecular  electric  and/or  magnetic  dipole  moments  in 
a  real  dielectric  and/or  ferrite  medium  [2].  These  macroscopic  dipole 
moments  modify  the  net  electric  and/or  magnetic  dipole  moment  per 
unit  volume,  and  thus  they  also  modify  the  effective  permittivity  and 
permeability  of  the  medium. 

An  interesting  and  potentially  useful  feature  of  artificial  media  is 
that  by  adjusting  the  size,  shape,  material  composition,  and  density 
of  the  scatterers,  one  may  be  able  to  engineer  or  design  a  medium 
which  has  desirable  permeability,  permittivity,  and  dispersion  charac¬ 
teristics.  Also,  the  theory  presented  here  may  be  useful  in  modeling 
conqxisite  materials  formed  by  placing  lossy  dipole  filaments  or 
periodic  weaves  of  long  lossy  filaments  in  a  host  medium.  The 
use  of  artificial  dielectrics  was  first  suggested  by  Kock  for  making 
microwave  lenses  [3].  Brown  [4]  showed  that  a  two-dimensional 
periodic  array  of  conducting  wires  is  an  artificial  dielectric  with 
index  of  reffaction  less  than  unity,  and  Bahl  and  Gupta  [5]  used  this 
idea  in  the  design  of  a  leaky  wave  antenna.  King,  Thiel,  and  Park 
inserted  small  pins  in  a  grounded  dielectric  substrate  to  synthesize  an 
artificial  dielectric  with  a  given  surface  reactance  [6].  Artificial  media 
have  also  been  used  in  the  design  of  microwave  lenses  [7}-{9],  by 
Rotman  to  model  plasmas  [10],  and  by  Sifavola  to  model  rain  and 


Manuscript  received  February  26,  1992;  revised  Sqiteinber  20.  1993.  This 
work  was  supported  by  die  Joint  Services  Electronics  IVogram  under  Contract 
N00014-78-C-0049  with  the  Ohio  State  University  Research  Foundation. 

The  authors  are  with  the  ElectroScience  Laboratory,  Ohio  State  University 
Columbus,  OH  43212  USA. 

IEEE  Log  Number  94013SS. 


hail  [11].  Hall  et  al  have  considered  the  use  of  artificial  dielectrics 
in  the  reduction  of  the  radar  cross  section  of  vehicles  [12]. 

Most  of  the  early  methods  of  analysis  for  artificial  media  were 
based  upon  static  approximations,  and  often  included  only  the  lower 
order  multipoles  [1,  Ch.  12].  For  diis  reason,  they  were  limited  to 
electrically  small  scatterers,  which  were  not  placed  too  close  to  one 
another.  More  recently,  quasistatic  solutions  have  been  presented  by 
Poulaiikas  [13]  and  Sihvola  [11].  By  contrast,  the  present  method 
is  a  full  wave  analysis  based  upon  a  method  of  moments  (MM) 
[14]  solution  of  an  exact  integral  equation  for  the  equivalent  currents 
representing  the  scatterers.  This  MM  solution  accounts  for  mutual 
coupling  between  the  scatterers,  is  ^iplicable  to  artificial  media 
comprised  of  scatteren  of  essentially  aibitrary  size,  shape,  and 
material  composition,  and  is  capable  of  determining  the  dispenion 
characteristics  of  the  artificial  medium.  As  opposed  to  many  other 
methods  [1,  Ch.  12],  it  is  not  necessary  to  know  the  electric  and 
magnetic  polarizability  of  the  scatteren,  since  this  is  inqilicitly 
determined  as  part  of  the  MM  solution.  In  this  and  most  other  amdyses 
of  artificial  media,  the  scatteren  are  assumed  to  be  identical,  and  to  be 
arranged  in  a  periodic  lattice.  This  simplification  allows  the  artificial 
media  to  be  analyzed  by  the  periodic  moment  method  (PMM)  [15], 
in  which  the  unknowns  in  the  MM  solution  are  limited  to  currents 
on  the  center  element 

Section  n  presents  the  general  integral  equation  and  PMM  analysis 
of  periodic  artificial  media,  and  includes  a  discussion  of  the  compu¬ 
tation  of  the  effective  permittivity  and  permeability  for  anisotropic 
artificial  media.  Numerical  results  for  a^cial  media  comprised  of 
a  3-D  array  of  short  dipoles  are  presented  in  Section  m,  and  are 
compared  with  static  solutions. 

n.  PMM  Analysis  of  AKnFiOAL  Media 
A  Overview  the  Problem 

As  illustrated  in  Fig.  1,  tiie  artificial  medium  is  modeled  as 
a  triple  infinite  periodic  array  of  identical  scattering  elements  in 
some  homogeneous  and  isotrc^ic  background  or  host  medium  with 
permittivity  and  penneability  denoted  {no,(o),  wavelength  Ao,  and 
wavenumber  ko  =  Although  here  no  assumption  is  made 

as  to  size  and  density  of  scatterers,  typically  tiiere  will  be  a  large 
number  of  electrically  small  scatterers  per  cubic  wavelengtii.  The 
spacing  between  the  elements  is  denoted  by  di ,  dj,  and  da  in  the  x,  y, 
and  z  directions,  respectively.  The  elements  in  the  triple  infinite  array 
of  scatterers  will  be  identified  by  the  integer  indices  Q  =  (/,m,n), 
(~o®  <  (].*»»>»»)  <  oo).  Element  Q  =  (0,0,0)  will  be  refeired  to 
as  the  center  element,  and  is  assumed  to  be  located  in  the  vicinity  of 
flic  origin.  AQ  =  {ldi±,md2y,ndix)  is  a  vector  firom  a  reference 
point  on  flie  center  element  to  flic  conesponding  point  on  element  Q. 
In  this  pqiCT,  all  fields  and  currents  are  assumed  to  be  time  harmonic, 
with  the  e'"''  time  dependence  siqipressed. 

The  basic  problem  is  to  determine  whether  or  not  the  periodic 
medium  of  Rg.  1  appears  as  an  artificial  medium,  and  ff  w,  to 
conqwte  the  effective  complex  permittivity  and  penneability  of  the 


0018-926X/94$(M.(X)  O  1994  TPRF 


728 


IEEE  TRANSACnONS  ON  ANTENNAS  AND  PROPAGAnC»I.  VCX.  42,  NO.  S.  HAY  1994 


Fig.  I.  A  triple  infinite  periodic  array  of  scatterers  tepre^ting  an  artificial 
dielectric. 


artificial  medium.  This  will  be  done  by  determining  whether  or  not 
a  plane  wave  is  an  eigenfunction  [16,  Ch.  4-1]  of  the  periodic 
medium  of  Fig.  1,  and  if  so,  what  is  the  eigenvalue  or  complex 
phase  constant  of  that  plane  wave.  By  an  eigenfunction  it  is  meant 
a  source-free  solution  of  Maxwell’s  equations  that  satisfies  all  of  the 
boundary  conditions.  To  find  the  eigenfunctions,  a  plane  wave  of  a 
known  frequency  w  is  assumed  to  be  propagating  through  the  periodic 
medium  in  a  known  direction.  In  a  source-free  homogeneous  medium 
(isotropic  or  anisotropic  [17,  Ch.  2.3]),  a  plane  wave  propagating  in 
the  known  u  (u  assumed  real)  direction  will  have  spatial  variation 
of  the  form 


C  9 


(I) 


where  k«  =  kcU  is  the  unknown  vector  phase  constant  of  the  plane 
wave  in  an  artificial  medium  of  wavenumber  and  r  =  rr  = 
+  3/y  +  **  is  a  vector  from  the  origin  to  the  point  (x,  y,  z). 
It  is  emphasized  that  at  this  point  in  the  solution  the  direction  of 
propagation  u  for  the  plane  wave  is  known;  however,  its  phase 
constant  kc  and  polarization  are  unknown.  Analogous  to  the  methods 
used  by  Itoh  and  Mittra  to  determine  the  effective  permittivity  of  a 
microstrip  transmission  line  [18],  if  a  value  of  kc  can  be  found  such 
that  this  plane  wave  is  an  eigenfunction  of  the  periodic  medium, 
then  the  periodic  medium  can  be  viewed  as  an  artificial  medium  with 
effective  dyadic  permeability  and  permittivity  denoted  (ii<,  («).  The 
evaluation  of  (fit,  it)  is  discussed  in  Section  II-C. 


B.  Integral  Equation  and  Moment  Method  Solution 

The  periodic  medium  of  Fig.  1  can  be  analyzed  using  the  periodic 
moment  method  (PMM),  which  has  been  extensively  applied  in  the 
analysis  of  periodic  arrays  or  frequency  selective  surfaces  [IS].  The 
solution  is  begun  by  using  the  equivalence  theorems  [2],  [16]  to 
replace  the  periodic  array  of  scatterers  by  the  host  medium  and  by 
equivalent  currents.  For  example,  the  volume  equivalence  theorem  [2, 
Ch.  7.7]  can  be  used  to  replace  dielectric  scatterers  by  electric  volume 
polarization  currents  located  in  the  volunte  occupied  by  the  scatterers. 
For  perfectly  conducting  scatterers,  the  surface  equivalence  theorem 
[2,  Ch.  7.8]  can  be  used  to  replace  the  scatterers  by  equivalent  electric 
surface  currents  on  the  surface  of  the  scatterers.  For  the  sake  of 
sinq)licity,  here  it  will  be  assumed  diat  the  equivalem  current  is 
electric  (as  opposed  to  magnetic),  and  thus  the  equivalent  electric 


current  on  the  elenoents  of  the  periodic  array  can  be  written  as 

(2) 

l,m,n 

where  is  the  equivalent  electric  current  on  element  Q,  and  the 
summation  is  over  all  array  elements  (-00  <  (l,m,n)  <  oo).  J** 
exists  over  the  region  (surface  at  volume)  of  element  Q.  Using 
the  periodicity  of  die  medium  and  of  the  plane  wave  field  of  (1), 
the  current  on  each  element  is  identical  except  for  an  amplitude  and 
phase  shift  corresponding  to  the  amplitude  and  phase  of  die  plane 
wave  at  the  reference  point  on  each  element  Using  (1),  the  complex 
tnagniimfe  of  the  current  on  element  Q,  relative  to  that  of  the  center 
element  is 


The  next  step  is  to  formulate  an  integral  equation  for  the  equivalent 
eigenfunction  current  J.  Although  there  are  many  possible  integral 
equ^ons  (surface,  volume,  electric  field,  magnetic  field,  etc.),  they 
can  always  be  written  in  the  form 

I(J)  =  0  over  R°  (4) 

where  £  is  a  linear  operator  and  is  the  region  over  which  the 
center  element  current  J°  is  defined.  Inserting  (2)  into  (4)  and  using 
the  linearity  of  L  yields 

L{3^)  =  0  over  R°.  (5) 

J,m,n 


In  order  to  solve  (5)  using  the  PMM,  the  current  on  die  entire 
array  is  expanded  as 

E  (6) 

{,m,n  j,m,n  jssl 


where  the  are  a  sequence  of  N  known  expansion  functions  on 
element  Q,  and  the  are  a  sequence  of  N  unknown  coefficients 
0  =  l,2,-”,iV).  The  expansion  frmcdons  jf  on  element  Q  are 
identical  to  those  on  the  center  element,  except  for  a  translation  of 
AQ. 

Inserting  (6)  into  (5),  and  using  the  linearity  of  L,  yields 

over  U®.  (7) 

/»m,n  7=1 

The  final  step  in  the  PMM  solution  requires  the  selection  of  N 
lineariy  independent  vector  weighting  functions  over  .R®,  denoted 
Wi,  i  =  l,2,”-,iV.  Taking  the  inner  product  of  (7)  with  each 
weighting  frinction  will  reduce  (7)  to  a  system  of  N  simultaneous 
linear  equations  in  the  AT  unknowns  Ij  (j  =  1, 2,  •  •  • ,  N),  which 
can  be  compactly  written  as  the  matrix  equation 


l^(k.)]/  =  0 


(8) 


where  \Z(kt)]  is  the  order  N  uripedance  matrix,  and  7  is  the  length 
N  current  vector  containing  the  N  coefficients  in  the  expansion  fm’ 
the  current  A  typical  element  of  the  impedance  matrix  is 

Zii(kt)  =  Y 

f,m,n 

=  C‘*(I(J?),W,)  =  (9) 

/,m,n 

where  Z,^  is  the  mutual  impedance  between  expansion  function  j  on 
element  Q  and  weighting  function  i  on  the  center  element,  and  the 
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brockets  indicate  an  appropriate  inner  product  Topically,  this  »n»w^ 
product  will  be  chosen  as 

f  I(jQ).w,dr  (10) 

Jro 

wdiere  die  integration  is  over  the  region  (line,  surface,  volume)  of  die 
center  element  for  which  w,  is  nonzero.  Note  that  the  of 

the  impedance  matrix  depend  upon  the  unknown  phase  constant  ik«, 
which  enters  through  the  of  (3).  The  impedance  sums  in  (9)  are 
typically  not  absolutely  convergent,  and  some  technique  to  accelerate 
convergence  is  usually  required  for  numerical  eompiitarinp^  []9]. 

The  homc^eous  matrix  equation  (8)  will  only  have  a  nontrivial 
solution  if 


l^(fc.)l  =  0, 


(11) 


that  is,  if  die  determinant  of  the  impedance  matrix  is  zero.  Equa¬ 
tion  (11)  must  be  solved,  usually  in  an  iterative  fashion,  for  the 
eigenvalues  k^,  which  result  in  a  zero  determinant.  If  a  root  of 
this  characteristic  equation  exists,  then  the  periodic  nwtiHtn  can 
be  viewed  as  an  artificial  medium,  with  effective  permeability  and 
permittivity,  as  described  in  the  next  section. 

To  find  the  eigenfunction  currents,  to  within  a  constant,  an  arbitrary 
(nonzero)  element  of  the  7  vector  is  set  to  unity,  and  (8)  is  reduced 
to  an  order  N  —  1  matrix  equation  for  the  remaining  —  1  elements 
in  7.  For  example,  setting  Is  =  1,  (8)  becomes 


[Z’{k.)]I'  = 


r  -Zis 
—Zis 

L->^n-i,n  J 


(12) 


where  the  [Z'\  matrix  holds  the  first  JV  —  1  rows  and  columns  of 
the  [Z]  matrix,  and  the  I'  vector  holds  the  first  AT  —  1  elements  of 
the  7  vector. 

The  eigenfunction  field  is  (voduced  by  die  superposition  of  die 
fields  of  the  currents  (— oo  <  (/,  m,n)  <  oo)  radiating  in  die 

homogeneous  host  medium  (f<o,  eo).  This  field  will  ^ipear  as  a  plane 
wave,  with  the  spatial  wiation  of  (1),  plus  some  fine  structure  or 
ripple  as  the  field  point  rqiproaches  one  of  the  scatterers. 


C.  Determination  of  the  Effective  Permeability  and  Permittivity 

Discussion  of  the  roots  and  polarization;  The  first  step  in  de¬ 
termining  the  equivalent  permeability  and  permittivity  of  the  artificial 
medium  is  to  determine  the  roots  from  (11).  For  a  given  direction 
of  propagation  through  an  anisotropic  artificial  nwtimn  composed 
of  arbitrary  scatterers,  diete  are  two  fundamental  roots  kt  to  (11) 
[20,  Sec.  4.25],  [21,  Sec.  14.2.2].  Each  root  cmresponds  to  a  distinct 
polarization  of  die  plane  wave  in  the  artificial  medium,  and  its 
corresponding  eigenfimction  currents  and  fields. 

For  special  geometries,  tbe  roots  may  be  rqieated  roots  or  degener¬ 
ate  to  the  host  medium  wavenumber  kt  —  ko.  For  example,  repeated 
roots  (two  roots  with  the  same  numerical  value  for  kt)  wiU  occur 
for  propagation  normal  to  a  symmetric  trite  cross  with  equal-length 
vemcal  and  horizontal  members.  One  root  corresponds  to  vertical 
polarization,  and  die  odier  root  corresponds  to  horizontal  polarization. 
If  the  scatterer  is  of  2rD  or  of  3-D  extent,  tbe  two  roots  will  be 
different  from  the  host  medium  wavenumber.  If  the  scatterer  is  of 
1-D  extent  (Le.,  a  linear  dipole),  dien  there  will  be  one  root  k.  ==  ko, 
corresponding  to  a  plane  wave  with  polarization  popendicular  to  the 
dipole,  and  one  root  kt  #  ko,  corresponding  to  polarization  parallel 
to  the  dipole. 

Note  that  in  developing  (11),  the  pcdarization  of  the  plane  wave 
was  not  qiecified.  Thus,  when  a  root  of  (1 1)  is  found,  the  polarization 
is  unknowtL  To  firul  the  ptdarization  one  can  first  compute  the 


correqxmding  eigenfunction  currents  using  (12),  and  then  die  rigen- 
functjon  electric  and  magnetic  fields.  The  polarization  of  the  electric 
field  is  taken  as  the  polarization  cmreqxinding  to  the  chosen  value 
of  k(.  Once  the  eigenfunction  fields  are  known,  the  characteristic 
impedance  of  the  artificial  medium,  denoted  tjt,  can  be  evaluated  as 
the  ratio  of  the  electric  to  magnetic  eigenfunction  fields  tangential 
to  the  assumed  direction  of  propagation.  If  the  host  medium  and  the 
scatterers  are  lossless,  then  kt  and  rjt  will  be  positive  real  numbers. 
However,  if  either  the  host  medium  or  the  scatterers  have  loss,  then 
kt  will  be  a  complex  number  in  the  fourth  quadrant,  and  will  be 
a  complex  number  in  tbe  sector  ^45°  of  die  positive  real  axis  [16, 
Sec.  2-3]. 

Since  the  scatterers  in  an  artificial  medium  ate  typically  electrically 
very  small,  one  is  usually  interested  in  die  lowest  order  eigenfunction 
modes  with  die  smallest  kt.  However,  it  is  irrqiortant  to  note  that 
larger  roots,  corresponding  to  higher  order  modes,  may  exist  In 
general,  we  know  of  no  a  priori  method  for  distinguishing  the  roots 
coiresponding  to  the  lowest  wder  modes  from  those  corresponding 
to  the  higher  order  modes.  Normally,  this  should  not  be  a  problem, 
since  for  the  electrically  small  scatterers  commonly  used  in  artificial 
media,  the  higher  order  modes  are  so  far  from  resonance. 

Anisotropic  artificial  media:  Cmsider  the  determination  of 
the  dyadic  effective  permeability  and  permittivity  {ptCt)  for  an 
anisotropic  artificial  medium.  It  is  mentioned  that  artificial  media  can 
be  chiral  or  even  bianisotropic  [17];  however,  (3)  would  no  longer  be 
valid,  and  this  is  beyond  the  scope  of  this  article.  It  is  assumed  that  for 
a  given  geometry,  tbe  two  roots  k,,  their  corresponding  eigenfunction 
currents  J**  on  the  center  element,  and  tlieir  eigenfunction  fields 
(E“,H®),  averaged  over  the  volume  of  die  center  cell,  have  all 
been  determined.  In  tbe  limit  as  ke  — *  ko,  the  eigenfunction  currents 
vanish  and  tbe  eigenfunction  fields  become  identical  to  a  plane  wave 
in  free  space.  Note  dial  in  anisotropic  media  (E'’,H°)  are  not  in 
general  ottbogonal  to  k«  [17,  Cb.  2.3]. 

The  dectric  and  magnetic  dipole  moment  per  unit  volume  in  tbe 
center  cell  can  now  be  computed  as  [1,  Cb.  12.5] 

where  71°  is  the  region  of  die  center  cell  of  volume  Ar  =  did^di, 
and  X*  and  x™  are  the  dimensionless  symmetric  [20,  Sec.  4.24], 
[21,  Sec.  14.1]  dyadic  effective  dectric  and  magnetic  susceptibilities, 
respectively,  for  the  artificial  medium.  Equation  (14)  shows  that  even 
perfecdy  conducting  or  pure  dielectric  scatterers  can  have  a  magnetic 
moment,  and  thus  an  effective  permeability  different  from  tbe  host 
medium.  Assuming  diat  die  usual  constitutive  relationships,  which  are 
valid  point  by  point  in  a  real  medium,  bold  in  an  average  sense  in 
an  artificial  medium,  the  average  dectric  and  magnetic  flux  densities 
in  tbe  center  cell  are  given  by 

D°  =  eoE°  +  P®  =  eoE®  +  eoX*  •  E°  =  e.  •  E®  (15) 

B°  =  ri-  M°)  =  ,4o(H°  4-  r  •  H°)  =  /i.  ■  H°.  (16) 

Bom  (15)  and  (16),  die  d)^c  effective  permittivity  and  permeability 
are  given  by 


?.  =  (l+X*)«o  (17) 

/i.  =  (I+x’")mo  (18) 

where  I  is  the  unit  dyad. 

Explichly  showing  the  dyadic  equation  (13)  rdating  the  effective 
dectric  susceptibility  to  die  average  dectric  field  and  die  dectric 
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dipole  moment  per  unit  vojume  in  the  center  cell. 
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Equation  19  i$  equivalent  to  three  equations  in  the  nine  components 
of  x',  and  is  the  result  of  one  of  the  two  roots  of  (11).  The  other 
root  will  produce  a  dyadic  equation,  similar  to  (19)  with  the  same 
x‘,  but  with  different  E°  and  P°.  The  two  dyadic  equations,  along 
with  the  condition  that  x‘  is  symmetric,  can  now  be  solved  fm  the 
nine  components  of  The  determination  of  fi,  is  parallel  to  that 
presented  for  ■  In  a  real  anisotn^ic  medium  the  elements  of  the 
dyadic  constitutive  parameters  are  independent  of  the  direction  of 
propagation  [20,  Ch.  4],  [21,  Sec.  14.1];  however,  this  may  not  hold 
in  an  artificial  medium. 

Uniaxial  artificial  media:  For  certain  sinqile  anisotropic  media, 
termed  uniaxial  media,  the  off  diagonal  components  of  x*  (and  x”*) 
are  negligible,  and  the  diagonal  components  are  related  to  E°  and 
P”  by 


Xii  “ 


E? 


<0 


i  =  z,y,  z. 


(20) 


Equation  20  can  be  used  to  find  the  tt  component  of  x‘  provided  Ef 
is  nonzero.  If  for  a  given  direction  of  propagation  and  root,  #  0 
for  i  =:  z,  y,  z,  then  (20)  can  be  us^  to  find  all  three  diagonal 
components  of  x*.  If  for  a  given  direction  of  propagation  and  root 
Ef  =  0,  then  if  =  0,  and  (20)  is  indeterminate.  In  this  case  x« 
can  be  determined  from  the  other  root,  or  by  a  different  direction  of 
propagation. 

Consider  propagation  along  of  one  of  the  three  principle  axes  in  a 
uniaxial  media.  The  two  roots  will  correspond  to  polarizations  in  the 
directions  of  the  two  axes  transverse  to  the  direction  of  propagation. 
For  this  special  case,  the  effective  permittivity  and  permeability  can 
be  determined  from  the  root  kt  and  the  characteristic  impedance 
r/e.  It  is  assumed  that  and  rjc  have  been  determined  for  a  given 
frequency,  polarization,  and  direction  of  propagation  along  one  of  the 
principle  axes.  Since  kt  and  are  related  to  and  tc  by 


kt  =  and 


(21) 


dien 


ke  ,  tfgkf 

ft  = -  and  lit  =  - — . 

urjt  ui 


U  fit  =  fio,  then  e,  is  given  by 


(22) 


ft  = 


kl 


w-'/io 

thus  eliminating  die  need  to  compute 


(23) 


m.  Numerical  Results 

This  section  will  present  numerical  results  for  artificial  media 
comprised  of  a  3-D  array  of  short  perfectly  conducting  z  directed 
dipoles  of  length  L  and  radius  a.  The  MM  solution  for  the  3-D 
dipole  array  is  based  upon  the  usual  thin-wire  formulation,  employing 
piecewise  sinusoidal  (PWS)  expansion  and  weighting  fimctions  [22], 
In  all  cases  the  background  medium  is  free  space  and  the  polarization 
of  the  electric  field  is  z.  For  this  case,  the  artificial  media  are  lossless 
artificial  dielectrics  (jit  =  /to)  with  teal  Ce  related  to  Jbc  by  (23).  The 


Fig.  2.  Relative  effective  permittivity  versus  lattice  size  for  a  3-D  array  of 
dipoles. 


Fig.  3.  Dispersion  of  the  effective  permittivity  for  a  3-D  array  of  dipoles 
(propagation  along  the  z  axis). 


artificial  media  ate  uniaxial,  with  ft  being  the  zz  component  of  ft. 
To  the  approximations  used,  the  zz  and  yy  components  ate  equal  to 
<0.  and  all  off  diagonal  components  are  zero. 

The  insert  in  Fig.  2  shows  a  3-D  array  of  short  dipoles  of  radius 
a  =  O.OOlAo  and  length  L  =  0.2Ao  in  a  lattice  with  spacing  di  =  dj 
and  is  =  di  -I-  L.  The  figure  shows  the  relative  effective  permittivity 
of  the  medium  as  a  function  of  lattice  size  di.  Curves  are  shown 
for  iV  =  1, 7,  and  15  PWS  basis  functions,  and  are  compared  to 
the  static  qtproxiroation  of  Collin  [1],  As  expected,  the  effective 
permittivity  increases  as  the  density  of  the  dipoles  increases.  Note 
that  for  relatively  large  lattice  spacing  the  current  on  the  short  dipoles 
is  nearly  triangular,  and  is  well  ^proximated  by  IV  =  1  sinusoidal 
basis  frmction  solutions.  However,  as  di/Ao  — »  0,  the  dipoles  are 
almost  touching.  This  causes  the  dipole  current  to  take  on  a  more 
*^at  top”  shape,  and  requires  mote  than  N  =  1  basis  fimctions  to 
obtain  a  converged  result 

The  final  set  of  data  will  illustrate  dispersion  in  artificial  dielectrics. 
The  insert  in  Fig.  3  shows  an  artificial  dielectric  comprised  of  dipoles 
of  length  L  =  0.6  cm,  radius  a  =  0.01  cm,  and  with  a  lattice  gracing 
di  —  d2  =  0.1  cm  and  da  =  0.7  cm.  For  a  wave  propagating  along 
the  X  axis,  the  =  1, 5,  and  11  term  MM  solutions  for  the  relative 
effective  permittivity  of  the  dipole  array  is  plotted  from  0.24  GHz  to 
24  GHz,  corresponding  to  dipole  lengths  from  0.0048Ao  to  0.48Ao. 
At  low  frequencies,  ffie  MM  solution  approaches  the  static  result 
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(dashed  line)  of  trt  «  2.5.  However,  as  the  frequency  increases, 
the  magnitude  of  the  scattered  field  from  each  dipole  increases,  thus 
causing  the  effective  permittivity  to  significantly  increase. 
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A  Partitioning  Technique  for  the  Finite  Element  Solution  of 
Electromagnetic  SMtterlng  flrom  Electrically  Large 
Dielectric  Cylinders 

Roben  Lee  and  Veera  Chupongstimun 

Absiroet — A  flnitr  eIcnieDt  parddoiiiiig  idieiiM  hat  been  developed  to 
reduce  the  computadonal  eoets  of  modeliiic  electricaU}'  larfe  geometries, 
lo  the  partidoninK  tebeme,  the  CTUoder  it  divided  into  many  aecdons.  The 
finite  ^mmi  method  it  applied  to  each  aeetion  independent  of  the  other 
■ections,  and  then  the  toludont  in  each  aecdon  are  coupled  throu^  the 
ate  of  the  taagendal  field  continuity  condldont  between  ^Jacent  aecdons. 
Since  the  coupling  matrix  it  dgnifleandy  nnaller  «>»■"  the  original  finite 
dement  matrix,  it  is  expected  that  both  the  CPU  fime  and  memory  cocts 
can  be  dgnificandy  reduced.  The  parddoning  eeheme  it  coupled  to  the 
hymoment  method  to  account  for  the  boundary  truncation.  Numerical 
r^ts  are  presented  to  demonstrate  the  cffideacy  and  accurac)  of  tbe 
method. 

I.  Introduction 

The  finite  element  method  (FBM)  is  being  used  extensively  to 
perform  electromagnetic  modeling,  because  it  is  efficient  in  mod¬ 
eling  geometries  that  are  both  penetrable  and  inhomogeneous.  The 
efficiency  is  largely  due  to  the  fact  that  the  resulting  matrix  is 
sparse  Banded  matrix  solvers  [1]  have  traditionally  been  used  to 
solve  the  FEM  stiffness  matrix.  However,  when  the  electrical  size 
of  the  geometry  of  interest  is  large,  the  number  of  unknowns 
required  to  obtain  an  accurate  solution  may  be  so  large  that  the 
use  of  banded  solvers  becomes  compuutionally  infeasible.  The 
development  of  direct  sparse  matrix  solvers  such  as  the  nested 
dissection  method  [2]  has  been  a  major  area  of  research.  These 
methods  are  capable  of  handling  multiple  excitations  with  only  a 
small  increase  in  compuution  cosu.  The  computational  complexity 
of  sparse  matrix  methods  is  usually  less  than  that  of  banded  matrix 
methods.  For  example,  the  computation  cost  of  the  nested  dissection 
method  is  proportional  to  .V’  *,  while  the  computation  cost  for  the 
banded  matrix  methods  is  proportional  to  .V’  for  two-dimensional 
geometries.  The  direct  sparse  matrix  methods  rely  heavily  upon  the 
use  of  matrix  partitioning  or  reordering  schemes  to  reduce  the  number 
of  nonzero  entries  in  the  factorization  of  the  resulting  matrixes. 

In  numerical  electromagnetics,  several  researchers  have  developed 
methods  that  have  charaaeristics  very  similar  to  those  used  in  direct 
sparse  matrix  methods.  However,  rather  than  formulating  it  from  a 
mathematical  viewpoint  of  reorganizing  the  matrix,  the  researchers 
considered  it  from  a  physical  viewpoint  where  the  geometry  is 
divided  into  subsections.  Within  each  section,  a  numerical  solution  is 
generated  independent  of  the  other  sections.  The  coupling  between 
the  sections  is  then  performed  in  either  an  iterative  manner  [3], 
(4]  or  by  employing  rigorous  and  more  computationally  expensive 
procedures  [S],  [6].  It  should  be  noted  that  all  these  techniques  have 
been  developed  for  integral  equation  methods.  In  this  p^>er,  a  parti¬ 
tioning  technique  is  presented  for  the  FEM  to  solve  the  problem  of 
electromagnetic  scanering  from  a  two-dimensional  cylinder.  Like  the 
integral  equation  methods,  the  geometry  is  divided  into  subsections 
in  which  a  set  of  finite  element  solutions  is  generated  for  each  section 
independent  of  the  other  sections.  Tbe  sections  are  then  coupled 
together  to  form  a  sparse  matrix  that  is  significantly  snudler  than 
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Hg.  1 .  Cross  section  of  a  representative  dielectric  cylinder. 
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Fig.  2.  A  representative  partitioning  of  the  dielectric  cylinder. 

the  original  FEM  stiffness  matrix.  This  partitioning  technique  can 
be  coupled  to  a  boundary  truncation  method  such  as  the  bymoment 
method  [7]  to  solve  the  scattering  problem.  During  this  research,  it 
was  discovered  that  another  group  was  employing  a  very  similar 
method  for  waveguiding  geomeuies  [8],  It  is  expected  that  both 
methods  will  produce  the  same  efficiency. 

n.  FORMLT,ATtON 

Ut  US  consider  an  infinitely  long  isotropic  dielectric  cylinder  whose 
shape  and  material  properties  are  invariant  along  the  c  axis,  as 
shown  in  Fig.  1 .  The  incident  wave  is  also  assumed  to  be  invariant 
in  c;  therefore  the  fields  around  the  cylinder  can  be  decoupled 
into  TM{Et,  H,.  Hy)  and  TE(H,.  E,.  Ey)  polarizations.  In  this 
p^rer,  the  TM  polarization  is  described.  With  some  slight  modifica¬ 
tions.  the  formulation  for  the  TE  polarization  can  also  be  obtained 
To  simplify  the  notation,  let  us  define  E  to  be  E- .  E  satisfies  the 
generalized  Helmholtz  equation,  which  is  given  by 

^  ’ (^^e)+-'VE=0  (1) 

where  e*  =  e  -  It  is  assumed  that  the  exptj-/)  time 

dependence  has  been  suppressed. 

Initially,  the  formulation  is  presented  to  solve  a  boundary  value 
problem  in  which  the  Neumann  boundary  condition  is  specified  on 
dS.ri  in  Fig.  1.  Then  the  procedure  for  coupling  the  partitioning 
method  to  the  bymoment  method  is  described  to  solve  the  scanering 
problem.  In  the  partitioning  method,  the  cylinder  is  divided  into 
M  sections  as  shown  in  Fig.  2,  where  dSi.-  ■  ■  -  ■  ■  .dSsi.] 

denote  the  interior  boundaries  that  are  created  from  the  partitioning 
operation.  The  segments  dSo  and  dS^i  are  artificial  constructs  that 
allow  us  to  represent  the  first  and  last  sections  as  four-sided  sections 

A.  FEM  Solution  in  a  Single  Section 

Let  us  consider  the  mth  section  of  the  partitioned  cylinder  as  shown 
in  Fig.  3.  The  solution  in  this  section  can  be  obtained  from  FEM  if 
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Fig.  3.  The  mth  tection  of  the  cylinder.  The  direction  of  the  Ungeotial 
magnetic  fields,  h„i  and  are  defined  as  shown. 

the  correct  boundary  conditions  are  specified.  However,  although  the 
Neumann  boundary  condition  on  d5,ii  is  assumed  to  be  known,  the 
boundary  conditions  on  dSm-i  and  dSm  are  not  known  a  priori.  Let 
us  define  h„L  and  fimR  to  be  the  tangential  magnetic  field  on  dSm-i 
and  <f5m.  respectively,  for  the  mth  element.  These  two  quantities  can 
also  be  wrinen  in  terms  of  a  derivative  of  the  electric  field  E  where 
their  orientation  is  shown  in  Fig.  3.  The  tangential  magnetic  fields 
can  be  written  in  terms  of  an  infinite  sum  of  linearly  independent 
functions  (♦,,  j  =  1,  2.  3,  •  ■  •)  multiplied  by  unknown  coefficients 
a;-)  and 

hmL(t)  =  ~  fe[0,dLl  (2) 

1=1 


1  4f 


±JL. 

H  an 


Rg,  4.  The  boundary  value  problem  in  the  mth  section  is  broken  down  into 
the  three  boundary  value  problems  shown. 


<€[o,d«)  (3) 

where  d/dn  represents  a  directional  derivative  whose  direction  is 
outward  normal  from  a  point  on  the  surface  of  the  mth  section.  The 
variable  f  represents  a  parametric  mapping  of  the  (x,  y)  coordinates 
along  dSm-i  to  the  line  defined  by  the  endpoints  t  =  0  and  f  =  di, 
where  di  is  the  length  of  dSm-i-  A  similar  parametric  mapping  is 
made  for  dSm,  where  dR  is  the  length  of  dSm. 

In  the  numerical  implementation  of  (2)  and  (3)  the  infinite  sums 
are  truncated  to  7m  z.  and  I„r  for  dSm-i  and  dSm,  respectively. 
Using  the  properties  of  linearity,  we  can  determine  the  electric  field 
Em  in  the  mth  section  by  superimposing  the  three  boundary  value 
problems  shown  in  Fig.  4,  where  Em,  Em,  and  Em  represent  the 
electric  field  solution  for  the  corresponding  boundary  value  problems 
in  the  figure.  Note  that  the  boundary  conditions  for  £„  and  Em  are 
given  by  (2)  and  (3);  therefore,  these  field  solutions  are  given  by 

£mir,  v)  =  ^  y)  (4) 

■  =1 

£«(X,1,)=  (5) 

1  =  1 

where  Af'm  is  the  finite  element  solution  for  the  boundary  value 
problem  in  Fig.  4  associated  with  £„  except  that  4'i(<)  replaces 
hmt(f)  as  the  boundary  condition.  Similarly,  aI^„  is  generated  from 
the  replacement  of  hmRit)  by  9,{t)  in  the  boundary  value  problem 
associated  with  £„.  The  total  field  in  the  mth  section  can  thus  be 
written  as 

Emix,  y)  =  E^  +  £  al^'^Atm  +  b[’”^A!‘,„.  (6) 

■=1  1=1 

Equation  (6)  is  true  for  all  sections  except  the  first  and  last  ones. 
For  the  first  section,  the  Neumann  boundary  condition  on  dSo  is 
known.  Thus,  the  second  term  in  (6)  can  be  directly  incorporated 
into  Em-  Similarly,  for  the  last  section,  the  third  term  in  (6)  can  be 
incorporated  into  E^. 


B.  Coupling  the  Solution  Between  Sections 

Once  the  solutions  have  been  generated  for  every  section  as 
described  above,  it  is  necessary  to  couple  them  together  in  order 
to  determine  the  coefficients  oj"’^  and  The  coupling  is  accom¬ 
plish^  by  enforcing  the  continuity  of  the  tangential  magnetic  and 
electric  fields  at  the  interface  between  adjacent  sections.  From  (2) 
and  (3),  the  continuity  of  the  tangential  magnetic  field  is  specified  by 

tmR  ^(m+l)t 

<=i  «=1 

m  =  1,  2.-  -,.\f- 1.  (7) 

The  negative  sign  on  the  right  hand  side  of  (7)  is  due  to  the 
orientation  of  the  tangential  magnetic  field  quantities  in  Fig.  3. 
Due  to  this  continuity  condition,  it  is  ^>propriate  to  assume  that 
ImR  =  7(m+))t  s  7m.  The  sum  in  (7)  can  be  equated  term  by 
term  to  obtain 

= -bj”’  m  =  l.  2,---,A7-l.  (8) 

The  use  of  (8)  in  (6)  eliminates  half  of  the  unknowns.  We  determine 
the  remaining  unknowns  by  enforcing  the  continuity  of  the  tangential 
electric  field  on  dSm, 

Em{x,  y)  =  E^-^-iix,  y)  {x,  y)  e  dSm-  (9) 

In  general,  it  is  not  possible  to  choose  the  coefficients  al"’’  and  bj'"' 
to  satisfy  (9)  exactly.  Instead,  (9)  is  enforced  in  a  weak  sense,  i.e.. 

/  E„{x,  y)u!,.m  dl  s  f  Em-t-iix.  yU  j.mdl  (10) 

JaSm  JiSm 

where  ity.m  (y  =  1,  2,  •  •  • ,  7m )  represents  a  set  of  linearly  indepen¬ 
dent  weighting  functions  on  dSm.  Substimting  (8)  into  (6)  and  men 
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Fig.  S.  The  block  tri-diagonal  coupling  matrix. 


(6)  into  (10),  we  obtain 


L 


•=i 


+  2-  0|  A,  („41) 


w 

J>is„ 


Vj.rr,  dl 


(11) 


This  equation  is  valid  for  tn  =  1  to  m  =  A/  -  1,  where  it  is 
•ssumed  that  fcj°'  s=  0  and  =  0.  The  matrix  that  results  from 
( 1 1 )  can  be  ordered  in  such  a  way  that  it  forms  a  block  tridiagonal 
matrix  as  shown  in  Fig.  5,  where  the  submatrices  and 

Cm  are  full  and  have  dimensions  of  Im  x  x  and 

Jm  X  /m+1 ,  respectively.  The  vectors  fm,  each  of  which  has  a  length 
/m .  represents  the  right  hand  side  of  (11).  The  solution  of  this  matrix 
equation  can  be  efficiently  computed  with  a  banded  matrix  solver  or 
a  matrix  solver  specially  designed  for  matrices  of  this  type,  such  as 
the  Ricatti  transform  [9]. 


C.  Boundary  Truncation 

To  solve  the  scattering  problem,  the  interior  FEM  problem  must  be 
coupled  to  the  exterior  free  space  region.  In  this  paper,  the  bymoment 
method  [7]  is  used  to  perform  the  coupling.  The  bymoment  method 
relies  upon  the  superposition  of  a  multiple  number  of  finite  element 
solutions,  which  are  generated  with  specified  Dirichlet  boundary 
conditions.  The  number  of  FEM  solutions  is  equal  to  the  number 
of  expansion  functions  used  to  ^rproximate  the  boundary  solution. 
The  bymoment  method  can  be  modified  such  that  the  multiple  finite 
element  solutions  are  generated  from  specified  Neumann  boundary 
conditions.  Thus,  the  partitioning  technique  can  be  used  to  generate 
each  of  the  finite  element  solutions.  The  actual  implemenution  of 
the  bymoment  method  can  significantly  affect  the  efficiency.  In  the 
current  scheme,  the  partitioning  technique  is  used  to  find  the  FEM 
solution  for  the  entire  computational  domain  due  to  each  Neumann 
boundary  condition.  This  scheme  is  the  simplest  to  implement,  but 
it  is  compuutionally  inefficient  because  the  finite  element  solutions 
associated  with  the  interior  boundaries,  A,^„  and  A{^„,  must  be 
recomputed  each  time  the  partitioning  technique  is  applied.  A  much 
more  efficient  scheme  would  be  to  generate  FEM  solutions  for  the 
bymoment  method  on  each  individual  section.  Then  the  coefficients 
associated  with  the  bymoment  method  on  each  section  can  be  placed 
in  the  coupling  matrix.  The  solution  of  this  coupling  matrix  provides 


a  complete  solution  to  the  scattering  problem.  We  are  currently 
implementing  and  testing  this  scheme. 

m.  Numerical  Conshcrahons 

To  determine  the  efficiency  of  the  partitioning  technique,  we 
perform  an  asymptotic  analysis  to  determine  the  number  of  floating 
point  operations  required  for  this  method  and  compare  the  results  to 
the  standard  method,  in  which  the  FEM  solutions  are  computed  with 
a  banded  matrix  solver  based  on  LU  decomposition  In  addition,  the 
memory  requirements  are  compared.  The  analysis  is  performed  for  an 
electrically  large  two-dimensional  square  grid  consisting  of  bilinear 
square  elements  with  N  nodes  (one  unknown  per  node)  in  the  grid. 
To  simplify  the  asymptotic  analysis,  we  assume  Neumann  boundary 
conditions  are  enforced.  For  the  standard  method,  the  number  of 
floating  point  operations  for  the  LU  decomposition  is  proportional  to 
A'®,  while  the  backsolve  is  proportional  to  A'*^’.  In  the  partitioning 
technique,  the  square  grid  is  partitioned  into  M  uniform  sections.  In 
this  analysis,  it  is  assumed  that  as  the  electrical  size  of  the  square 
increases,  the  number  of  sections  M  iiKreases  in  such  a  way  as  to 
maintain  a  constant  width  for  each  section.  The  LU  decomposition 
time  associated  with  each  section  it  greatly  reduced  because  the 
bandwidth  of  the  FEM  matrix  remains  constant.  Consequently,  the 
sum  of  the  decomposition  times  for  all  the  sections  is  proportional  to 
A*.  The  partitioning  method  requires  performing  numerous  backsolve 
operations  on  each  section.  The  number  of  floating  point  operations 
required  to  perform  the  backsolves  for  all  the  sections  is  proportional 
to  A’®''*.  The  solution  of  the  final  coupling  matrix  shown  in  Fig.  5 
requires  an  operation  count  that  is  proportional  to  Af.V®''®.  Since 
Af  is  proportional  to  A’’^®,  the  complexity  of  the  partitioning 
method  is  equal  to  that  of  the  standard  method,  but  Af  is  usually 
chosen  to  have  values  of  A'’^®/20  to  A’’''®/40.  Thus,  the  operation 
count  can  be  reduced  by  a  large  multiplicative  factor.  However,  in 
practice,  the  cylinders  must  be  very  large— over  50  wavelengths 
in  both  directions — for  the  computation  time  associated  with  the 
coupling  matrix  to  dominate  the  overall  computation  time.  For  smaller 
cylinders,  we  expect  the  backsolve  time  to  be  the  dominant  factor. 

In  addition  to  improving  the  computational  efficiency,  the  parti¬ 
tioning  technique  can  substantially  reduce  the  memory  storage  from 
the  standard  jnethod.  The  memory  requirements  of  a  banded  matrix 
algorithm  is  A’®''®.  In  the  partitioning  technique,  each  section  can  be 
handled  in  a  sequential  manner;  therefore,  the  same  memory  locations 
can  be  used  for  each  section,  resulting  in  a  memory  requirement  of 
A’®^®/A/®.  The  extra  factor  of  Af  in  the  denominator  is  due  to  the 
fact  that  the  half-bandwidth  of  the  FEM  matrix  for  each  section  has 
been  reduced  by  a  factor  of  Af .  For  the  coupling  matrix  in  Fig.  5,  the 
memory  requirement  is  proportional  to  A'Af ,  where  the  exact  value 
is  dependent  upon  the  type  of  basis  functions  used  in  (2)  and  (3).  For 
the  problems  that  we  consider  in  this  p^>er,  a  factor  of  ten  reduction 
in  the  memory  is  typical. 

There  are  several  other  possible  advantages  of  the  partitioning 
technique.  For  problems  in  which  the  material  properties  change 
abruptly,  the  FEM  grid  must  transition  from  one  grid  density  to 
another  at  the  material  interface  in  order  to  maintain  a  constant  nodal 
density  per  wavelength  of  the  material.  This  transition  is  not  easy 
to  achieve  for  traditional  finite  element  techniques  However,  in  the 
partitioning  method  each  material  can  be  placed  in  a  separate  section 
Since  the  grid  for  each  section  is  independent  of  the  other  sections, 
there  is  no  need  to  provide  a  transition  between  the  sections 

In  partitioning  the  geometry  into  sections,  there  are  many  instances 
where  two  or  more  of  the  sections  are  geometrically  identical  In  these 
instances,  the  FEM  solution  need  only  be  computed  once  for  each 
set  of  identical  sections,  thereby  further  reducing  the  computational 
costs.  A  further  advantage  of  the  partitioning  technique  is  that  it  is 
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Ciisily  parallelized  for  multiprocessor  computers,  because  each  section 
can  be  assigned  to  a  processor  without  any  need  to  communicate 
between  the  processors  until  the  final  evaluation  of  the  coupling 
matrix.  The  final  coupling  matrix,  which  is  block  tiidiagonal  is  also 
highly  parallelizable  [lOJ. 


IV.  Numerical  Results 

To  demonstrate  the  efficiency  of  the  partitioning  technique,  we 
consider  two  boundary  value  problems  and  two  associated  scattering 
problems  The  excitation  for  all  the  geometries  is  assumed  to  be  a 
TM  polarized  plane  wave  of  the  form 

£;""  =  foe"-’*'”'”*'-’' (12) 


The  geometry  for  the  boundary  value  problems  consists  of  a  rect¬ 
angular  region  of  free  space.  This  geometry  was  chosen  because 
the  Neumann  boundary  condition  for  a  free  space  region  is  well 
known.  For  the  scattering  problems,  rectangular  dielectric  cylinders 
are  considered.  The  parameter  of  interest  is  the  echo  width  per 
wavelength  L,/X  [11],  and  it  is  given  by 


Ao 


lim 

p— 3C 


Ao 


(13) 


where  £*''  is  the  scanered  electric  field  and  Ao  is  the  free  space 
wavelength.  The  echo  width  is  computed  as  a  function  of  6  where  6 
is  the  angle  defined  from  the  positive  x  axis. 

A  finite  element  code  based  on  the  eight-node  quadrilateral  [1] 
element  has  been  written  to  test  this  method.  A  nodal  density  of 
approximately  20  nodes  per  wave-length  has  been  maintained  for  the 
numerical  computations.  This  choice  of  nodal  density  is  expected  to 
produce  accurate  results  (12).  The  geometries  of  interest  are  electri¬ 
cally  large;  therefore,  the  FEM  solutions  may  be  subject  to  errors 
produced  by  the  interior  resonances  of  the  compuUtional  domain 
[13].  To  remove  these  errors,  a  small  loss  has  been  introduced  in  the 
computation  domain.  Thus,  we  set  =  fr  —  )  10"*,  where  e,  is  the 
relative  permittivity  given  by  the  geometry.  The  computation  times 
shown  throughout  this  section  are  for  the  Cray  YMP  supercomputer. 
The  basis  and  weighting  functions  '1',  and  u'j.„  may  be  either  entire 
domain  functions  such  as  sinusoidal  functions  with  support  over  dSm 
(Fig.  1)  or  else  subdomain  functions  such  as  triangle  functions  with 
support  over  a  subsection  of  5m .  For  the  numerical  results  presented 
in  this  paper,  triangle  functions  with  support  over  two  elements  (5 
nodes)  are  used. 

The  first  example  is  the  problem  of  a  plane  wave  propagating 
through  a  2Ao  x  24Ao  region  of  free  space.  A  finite  element  mesh 
of  16141  nodes  (5200  elements)  is  used  to  discretize  the  entire 
compuution  domain.  In  Table  I,  a  comparison  is  made  between  the 
standard  method  and  the  partitioning  technique  in  which  the  domain 
has  been  divided  into  either  six  (2Ao  x  4Ao)  or  12  (2Ao  x  2Ao) 
sections.  Since  the  geometry  of  each  section  is  the  same,  only  one 
section  is  considered  in  the  FEM  computation.  From  the  table,  it 
is  evident  that  the  reduction  in  CPU  time  is  undramatic,  especially 
given  the  fact  that  the  FEM  solution  is  only  computed  in  one  section. 
The  improvements  that  we  expect  from  the  asymptotic  analysis  are 
not  present  in  this  example  because  the  bandwidth  of  the  FEM 
matrix  for  each  section  is  the  same  as  the  bandwidth  for  the  entire 
computation  domain.  Thus,  a  reduction  in  CPU  time  cannot  be 
expected  for  elongated  geometries  that  are  electrically  large  only 
along  one  direction.  The  reduction  in  memory  requirements,  on 
the  other  hand,  is  significant  (a  factor  of  more  than  eight  for  the 
12-section  case). 

The  second  example  is  a  Kattering  problem  based  on  the  first 
example.  The  scatterer  is  a  lAo  x  12Ao  ^electric  cylinder  (fr  —  4). 
The  mesh  is  the  same  as  the  one  used  in  the  first  example.  Because  the 


TABLE  I 

A  COMTAIUSON  OF  CPU  TIME  AND  MEMORY  REQUREMEKT  BETWEEN  THE 
Standard  Method  and  the  PARrmoNiNC  Technique  for 
four  Difperent  Examfles 


Example 

Number 

Solution  Method 

CPU  Time 
(Seconds) 

Memory 

(Megawords) 

1 

Standard 

16.2 

2.5 

Partitioning  (6  tections) 

9.2 

0.6 

Partitioning  (12  tectiont) 

7.4 

0.3 

2 

Standard 

156.5 

6.7 

Partitioning  (6  sections) 

221.7 

3.1 

Partitioning  (12  tections) 

286 

2.8 

3 

Standard 

388.5 

26.1 

Partitioning  (5  sections) 

34.5 

1.3 

Partitioning  (10  sections) 

27.6 

1.1 

4 

Standtud 

608.0 

28.3 

Partitioning  (5  tections) 

287.5 

3.0 

Partitioning  (10  sections) 

274.3 

2.8 

i 
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Fig.  6.  Plot  of  the  echo  width  as  i  function  of  $  for  lAo  x  I2Ao  cylinder 
with  C>,  =:  90°.  (r  =  4.  The  solution  from  the  standard  method  is  compared 
against  the  solution  from  the  partitioning  technique. 


geometry  is  electrically  large,  226  expansion  functions  are  required 
for  accurate  implementation  of  the  bymoment  method,  which  results 
in  a  dramatic  increase  in  coiiqiutation  costs.  The  echo  width  solution 
is  shown  in  Fig.  6  for  an  incidence  angle  of  ,90®.  The  results  from 
the  partitioning  technique  and  the  standard  method  are  vinually 
identical,  demonstrating  that  the  partitioning  has  very  little  efTeci  on 
the  solution.  Because  the  bandwidth  of  the  FEM  matrix  is  small,  the 
majority  of  the  CPU  time  is  spent  on  the  boundary  truncation.  From 
Table  I,  we  see  that  the  partitioning  technique  is  less  efficient  than  the 
standard  FEM  solution.  This  is  due  to  the  inefficient  manner  in  which 
the  bymoment  method  is  coupled  to  the  partitioning  technique  The 
efficiency  is  expected  to  improve  significantly  with  a  better  scheme. 
The  memory  savings  has  also  been  reduced  compared  to  the  first 
example  due  to  the  memory  requirements  of  the  bymoment  method 
For  the  third  example,  we  consider  a  plane  wave  propagating 
through  a  9Ao  x  9Ao  region  of  free  space.  Because  the  region  is  a 
square,  it  is  expected  to  exhibit  the  computational  savings  described  in 
the  asymptotic  analysis.  The  computation  domain  is  discretized  with  a 
grid  consisting  of  29  105  nodes  (9580  elements).  Fur  the  partitioning 
technique,  the  domain  is  divided  into  both  five  (9Ao  x  l.SAo)  and 
ten  (9Ao  x  0.9 Ao)  sections.  In  Table  I  the  computation  times  are 
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Fig  7.  Plot  of  the  echo  width  as  a  function  of  6  for  4.5Ao  x  4.5Aa  cylinder 
with  s  45°.  (r  =  4.  The  solution  from  the  standard  method  is  compared 
against  the  solution  from  the  partitioning  technique. 


tabulated  for  all  three  cases.  The  reduction  in  CPU  time  is  very 
significant  (more  than  a  factor  of  14  for  the  ten-section  case).  Again, 
we  took  advanuge  of  the  fact  that  the  FEM  solution  only  needs  to  be 
computed  in  one  section  since  all  the  sections  are  identical.  However, 
if  the  compuution  is  performed  for  all  the  sections,  the  CPU  time 
increases  from  34.5  s  (for  the  five-section  case)  to  ^rproximately  65 
s.  Thus,  for  square  dielectric  cylinders  with  nonidentical  sections,  the 
compuutional  savings  is  still  very  large.  Furthermore,  the  memory 
requirements  have  been  reduced  from  26  Mwords  for  the  standard 
method  to  1.1  Mwords  for  the  partitioning  method. 

In  the  fourth  and  final  example,  the  plane  wave  scattering  from 
a  square  dielectric  (<r  =  4)  cylinder  with  sides  of  length  4.5Ao  is 
considered.  The  FEM  grid  is  the  same  as  the  one  used  in  the  third 
example,  and  the  bymoment  boundary  condition  is  implemented  with 
180  expansion  functions.  The  echo  width  solution  is  shown  in  Fig.  7 
for  di  =  45°.  The  two  methods  produce  virtually  identical  resulu.  A 
comparison  of  the  memory  and  CPU  requirements  for  this  example 
is  shown  in  Table  I.  For  the  partitioning  technique,  the  compulation 
domain  is  divided  into  both  five  and  ten  sections.  The  reduction  in 
CPU  time  is  not  as  dramatic  as  in  the  third  example  because  of  the 
inefficient  manner  in  which  the  bymoment  method  is  coupled  to  the 
partitioning  technique.  However,  the  savings  are  still  very  significant 
in  both  computation  time  and  memory  costs. 


REfERENCES 

(1]  ).  E.  Akin,  Application  anJ  iii^lemeiuaiwn  of  Finite  Element  Methods. 
San  Diego;  Academic,  1982. 

(2]  A.  Ccorge  and  ).  W.  Liu,  Computer  Solution  cf  Large  Sparse  Positive 
Definite  Systems.  Englewood  Cliffs.  NJ:  Prentice-Hall.  1981. 

(3]  K.  R.  Umaahankar  and  A.  Taflove,  *7<lumerical  analysis  of  electromag¬ 
netic  icanering  from  electrically  large  objects  using  spatial  decomposi- 
two  lechniques,”  /£££  Trans.  Antennas  Propagoi. .  vol.  40.  pp.  867-877 
Aug  1992. 

(4]  C.  M.  Butler,  *Diakoptic  theory  and  the  moment  method."  IEEE 
Aruennas  Propagoi.  Soc.  Irtt.  Sym.  Dig.,  vol.  l.pp.  72-75.  May  1990. 

{5)  W.  C.  Chew  and  C.  C.  Lai,  “NEPAL- An  S'  °  algorithm  for  solving 
volume  integral  equations,"  IEEE  Antennas  Propagoi.  Soc.  Int  Sym. 
Dig.,  vol.  1.  pp.  184-187.  July  1992. 

(6)  T.  M.  Wang  aixi  H.  Ling,  *^ectromagnetic  icanering  from  three- 
dimensional  cavities  via  a  connection  scheme."  IEEE  Trans  Antennas 
Propagat.,  vol.  39.  pp.  1505-1513,  Oct.  1991, 

(7)  A.  C.  Cangellaris  a^  R.  Lee,  The  byirwment  method  for  two- 
dimensional  electromagnetic  tcattering."  IEEE  Trans.  Antennas  Prop¬ 
agat.,  vol.  38.  pp.  1429-1437,  Sept  1990. 

(8)  C.  T.  Spring  and  A.  C.  Cangellaris,  "A  partitioning  approach  to  the 
clectrornagnetic  characterization  of  large  two  and  three  dimensional 
structures,"  1992  URSI  Radio  Science  Meeting  Digest,  p.  230. 

{9)  K.  K.  Mei,  “Unimotnent  method  of  solving  antenna  and  scattering 
problems,"  IEEE  Trans.  Antennas  Propagat.,  vol.  22,  pp.  760-766.  Nov. 
1974. 

(10)  M.  Berry  and  A.  Sameh,  “Multiprocessor  schemes  for  solving  block 
Iridiagonal  linear  systems,"  Int.  J.  Supercomputing  Appl,  vol.  2,  pp 
37-57.  1988. 

(11)  9..  P.H»st\sti^on,Tmte-Honnonic  Electromagnetic  Fields  New  York: 
McGraw-Hill.  1961, 

(12)  R.  Lee  and  A.  C.  Cangellaris.  “A  study  of  discretization  error  in  the 
finite  element  approximation  of  wave  solutions."  IEEE  Trans.  Antennas 
Propagat.,  vol.  40.  pp.  542-549,  May  1992. 

(13)  L.  W.  Pearson,  A.  F.  Peterson,  L.  J.  Bahrmasel,  and  R  A  Whitaker, 
“Inward-looking  and  outward-looking  formulations  for  scattering  from 
penetrable  objects,"  IEEE  Trans.  Antennas  Propagat.,  vol.  40.  pp. 
714-720,  June  1992. 


V.  Summary 

A  new  method  has  been  presented  to  decrease  the  solution  time  for 
the  finite  element  modeling  of  electromagnetic  scattering  problems. 
In  this  method,  the  geometry  is  partitioned  into  smaller  sections. 
Within  each  section,  a  set  of  finite  element  solutions  is  generated 
that  are  independent  of  the  other  sections.  The  sections  are  then 
coupled  by  the  use  of  the  field  continuity  conditions.  The  resulting 
matrix  equation  is  block  tridiagonal  and  significantly  smaller  than  the 
original  FEM  matrix  equation.  An  asymptotic  analysis  was  provided 
to  demonstrate  the  potential  efficiency  of  this  method.  The  method 
was  coupled  to  the  l>ymoment  method  in  order  to  solve  the  scattering 
problem.  Finally,  numerical  resulu  were  presented  to  demonstrate 
both  iu  accuracy  and  efficiency  for  some  canonical  geometries. 
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Abstract — The  riectromagnetic  modeling  of  engine  cavities  is  a 
ver>  difficult  task  because  the  electrical  size  of  the  cavity  is  veiy 
large,  while  the  engine  termination  is  geometrically  complex. 
High-frequency  techniques  can  adequately  model  the  cavity,  but 
perform  poorly  when  applied  to  the  termination.  Low-frequency 
techniques  are  currently  infeasible  for  such  large  geometries 
because  of  the  large  memory  and  computation  time  require¬ 
ments.  In  this  paper,  we  present  a  hybrid  method  which  com¬ 
bines  the  roost  attractive  features  of  the  low-  and  high-frequency 
techniques.  The  finite-difference  time-domain  (FDTD)  method  is 
applied  to  the  small  region  surrounding  the  termination.  The 
remainder  of  the  cavity  is  modeled  with  ray  methods.  To  validate 
this  method,  we  consider  two-dimensional  cavities  with  complex 
terminations.  Our  results  are  compared  against  those  found 
from  a  hybrid  combination  of  the  modal  method  and  the  method 
of  moments. 


I,  Introduction 

The  problem  of  electromagnetic  scattering  from 
an  open-ended  cavity  structure  is  an  important  prob¬ 
lem  in  the  area  of  radar  scattering.  Among  other  things,  it 
models  the  radar  scattering  from  a  jet  engine  inlet.  Tradi¬ 
tionally,  high-frequency  methods  have  been  applied  to 
solve  such  problems.  In  recent  years,  two  ray-based  meth¬ 
ods  have  been  widely  used  to  model  the  scattering  proper¬ 
ties  of  electrically  large  cavities  of  arbitrary  shape.  These 
are  the  “shooting  and  bouncing  ray”  (SBR)  method  [1] 
and  the  “generalized  ray  expansion”  (GRE)  method  [21 
However,  these  two  methods  have  only  been  shown  to 
work  for  cavities  with  planar  or  very  simple  terminations. 
The  terminations  of  practical  interest  are  often  geometri¬ 
cally  complex  (such  as  the  fan  blades  in  a  jet  engine),  and 
the  ray  ihethods  are  not  expiected  to  be  able  to  properly 
account  for  the  diffraction  effects  of  the  termination.  In 
addition,  the  termination  may  contain  materials  which  are 
not  perfectly  conducting.  Thus,  a  solution  based  purely  on 
ray  methods  may  be  drastically  different  from  the  correct 
solution.  Low-frequency  techniques  such  as  the  moment 
method  (MM)  [3],  the  finite-element  method  (FEM)  [4], 
and  the  finite-difference  time-domain  method  (FDTD)  [5], 
[6]  can  accurately  model  the  complex  termination  as  well 
as  the  rest  of  the  cavity,  but  unfortunately,  the  computa- 
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tion  costs  to  model  such  a  large  geometry  exceed  even  the 
most  powerful  supercomputer. 

In  this  paper,  we  present  a  hybrid  method  which  com¬ 
bines  both  high-  and  low-frequency  methods  to  overcome 
the  shortcomings  of  either  methods.  A  ray  method  is  used 
to  track  the  fields  from  the  mouth  of  the  cavity  up  to  an 
arbitrarily  defined  planar  surface  close  to  the  termination. 
These  fields  act  as  the  excitation  for  the  termination 
region,  which  is  modeled  by  a  low-frequency  method.  In 
our  case,  we  have  chosen  the  FDTD  method  because  of 
its  computational  efficiency  and  its  capability  to  obtain 
multiple-frequency  data  from  a  time  history  solution.  It  is 
expected  that  this  hybrid  method  can  provide  a  solution 
that  is  almost  as  accurate  as  if  the  solution  had  been 
obtained  from  low-frequency  techniques  but  at  a  very 
small  fraction  of  the  computation  costs. 

The  paper  is  divided  into  three  sections  after  this  intro¬ 
duction.  Section  II  describes  the  details  of  the  hybrid 
method  including  a  description  of  the  coupling  between 
the  ray  method  and  the  FDTD  method.  Although  the 
formulation  of  the  hybrid  method  is  for  the  general  case 
of  a  three-dimensional  (3D)  cavity  with  an  arbitrary  termi¬ 
nation,  the  validation  of  this  method  is  performed  in  two 
dimensions  (2D)  in  order  for  us  to  more  easily  analyze  the 
accuracy  and  robustness  of  this  method  before  implemen¬ 
tation  in  3D.  In  Section  III,  results  are  presented  for  the 
2D  geometry  of  a  parallel  plate  waveguide  containing 
perfectly  conducting  terminations.  These  results  are  com¬ 
pared  to  a  hybrid  modal-MM  method  to  demonstrate  the 
accuracy  of  the  hybrid  ray-FDTD  technique.  Finally,  a 
brief  summary  of  this  paper  is  provided  in  the  last  section. 

II.  The  Hybrid  Ray-FDTD  Method 

In  this  section,  a  brief  discussion  of  the  GRE  and  the 
FDTD  methods  is  given,  followed  by  a  description  of  how 
the  coupling  between  the  frequency-domain  GRE  method 
and  the  FDTD  method  is  accomplished.  The  use  of  ab¬ 
sorbing  boundary  conditions  and  the  computation  of  the 
cavity  scattered  field  are  then  considered. 

A.  The  CRE  Method 

SBR  and  GRE  are  the  two  ray  methods  that  have  been 
considered  for  use  in  combination  with  FDTD  for  the 
hybrid  method.  Many  papers  have  already  been  published 
on  the  two  ray  methods  and  their  applications  in  cavities. 
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comer  reflectors,  and  antennas  [7]-[12].  For  the  hybrid 
method,  we  have  chosen  GRE  over  SBR  after  due  consid¬ 
eration  of  their  accuracy  and  efficiency.  As  in  any  hybrid 
method,  the  accuracy  of  the  solution  at  the  interface 
between  the  methods  is  an  important  factor  in  the  overall 
accuracy  of  the  hybrid  method.  It  has  been  shown  in  [13] 
that  the  ray  fields  predicted  by  GRE  inside  a  semi  infinite 
parallel-plate  waveguide  cavity  is  more  accurate  than  those 
obtained  by  SBR.  The  difference  in  the  accuracy  of  these 
ray  methods  is  attributed  to  the  difference  in  the  coupling 
of  the  (externally  applied)  incident  fields  into  the  wave¬ 
guide  cavity.  GRE  intrinsically  includes  the  incident  fields 
diffracted  into  the  cavity  by  the  open  end  through  the 
integration  of  the  equivalent  currents  over  the  aperture  of 
the  cavity.  SBR,  on  the  other  hand,  does  not  include  the 
edge-diffracted  fields.  Depending  on  the  electrical  size  of 
the  aperture  and  the  length  of  the  cavity,  the  edge-dif¬ 
fracted  fields  may  eventually  affect  the  geometrical  optics 
fields  intercepted  by  the  aperture.  In  terms  of  efficiency, 
SBR  generally  requires  less  rays,  which  also  means  less 
ray  tracing  to  perform,  than  in  GRE  for  a  single  incident 
angle.  However,  a  new  set  of  rays  has  to  be  traced  in  SBR, 
including  computing  their  divergence  factors  and  reflec¬ 
tion  coefficients,  for  every  new  incident  angle,  whereas 
the  same  set  of  rays  is  used  in  GRE  regardless  of  the 
incident  angle  of  the  field  (the  latter  will  become  evident 
later).  So  GRE  is  better  suited  for  scattering  problems 
where  the  scattered  field  is  required  over  a  range  of 
angles.  Therefore,  for  accuracy  and  efficiency,  we  have 
chosen  GRE  over  SBR  for  the  hybrid  method. 

Because  the  detailed  formulation  of  the  GRE  method 
can  be  found  in  [2],  only  the  working  principles  of  the 
method  will  be  described.  In  this  method,  the  incident 
field  at  the  aperture  of  the  cavity  is  replaced  by  equivalent 
surface  currents  via  Kirchhoffs  approximation.  These 
equivalent  currents  radiate  the  desired  fields  into  the 
cavity  interior  as  shown  in  Fig.  1.  For  large  cavities,  the 
aperture  is  divided  into  a  number  of  smaller  subapertures. 
A  cone  of  ray  tubes  is  then  launched  into  the  cavity  from 
the  phase  center  of  each  subaperture.  The  ray  tubes  are 
tracked  using  the  central  rays  of  the  respective  ray  tubes. 
The  central  rays  are  tracked  as  ordinary  geometrical  op¬ 
tics  rays.  The  exiting  ray  tubes  contribute  to  the  scattered 
field  of  the  cavity.  The  scattered  field  is  computed  by 
physical  optics  approximation,  taking  into  account  the 
wavefront  curvature,  the  size,  and  the  shape  of  the  ray 
tube  associated  with  each  exiting  ray  [2],  [14].  A  more 
efficient  method  for  obtaining  the  scattered  field  will  be 
presented  in  a  later  section. 

In  GRE,  each  ray  tube  is  weighted  by  the  far-field 
approximation  of  the  radiation  pattern  of  the  subaperture 
(from  which  that  ray  tube  has  originated)  with  the  cavity 
walls  absent.  Specifically,  the  electric  field  along  the  pth 
ray  of  the  /th  subaperture  prior  to  any  reflection  is  given 
by 


£p,(r)  =  C,(v) 


exp  ( -jkotp,) 

^Pi 


(1) 


Fig.  1.  The  generalized  ray  expansion  method. 


where  is  the  propagation  constant  of  the  medium  in 
which  the  cavity  is  embedded,  r^,  is  the  position  vector  of 
a  point  along  the  pth  ray  with  respect  to  the  phase  center 
O,  of  the  /th  subaperture.  C,(rpi)  is  the  far-field  vector 
radiation  pattern  of  the  electric  field  evaluated  in  the 
direction  in  the  absence  of  the  cavity  walls.  It  is  given 
by 

c,(V)-^//Jvxvxy«(r;) 

+  Yorpi  X  Af/Hr;)]  exp  ■  r])  ds'„  (2) 

where  rj  is  the  vector  from  O,  to  the  equivalent  sources 
on  the  subaperture  S,.  Equation  (2)  can  be  evaluated 
numerically  or  in  closed  form  depending  on  the  field 
excitation  and  the  shape  of  the  subaperture  [9]. 

As  each  ray  tube  undergoes  multiple  reflections  within 
the  cavity,  its  field  amplitude  is  modified  by  divergence 
factors  and  reflection  coefficients.  In  particular,  the  elec¬ 
tric  field  Eir^)  of  the  ray  immediately  after  reflection 
from  the  point  r,  is  found  iteratively  from  the  field 
£(r,_,)  immediately  after  the  previous  reflection  point 
Tj.,  through  the  relation 

£(r,)  =  r-£(r,.,)(DF),-ie->*»%  (3) 

where  s  -  |r,  -  F  is  the  planar  dyadic  reflection 
coefficient  evaluated  at  and  (DF),_,  is  the  divergence 
factor  governing  the  spreading  of  the  ray  tube  associated 
with  the  ray  after  reflection  from  the  point  at  r,_,.  The 
divergence  factor  is  usually  evaluated  via  the  0-matrix 
formulation  of  Deschamps  [IS]. 

From  (l)-(3),  it  is  clear  that  only  one  set  of  rays  needs 
to  be  traced  for  all  the  incident  angles  because  the  ray 
divergence  factors  and  reflection  coefficients  are  functions 
of  the  ray  paths  and  not  the  incident  angle;  only  the  initial 
amplitude  of  the  ray  field  is  a  function  of  the  incident 
angle. 

The  most  serious  limitation  in  ray  methods  (like  the 
GRE  and  SBR  methods)  is  their  inability  to  account  for 
fields  arising  from  discontinuities  like  edges  and  tips  in¬ 
side  the  cavities.  As  such,  ray  methods  are  not  expected  to 
produce  reasonable  results  for  cavities  with  complex  ter¬ 
minations. 

B.  The  FDTD  Method 

FDTD  is  a  direct  solution  of  Maxwell’s  time-dependent 
curl  equations  [5],  [6].  It  applies  second-order  accurate 
central-difference  approximations  for  both  the  space  and 
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time  derivatives  of  the  electric  and  magnetic  fields  directly 
to  the  differential  operators  of  the  curl  equations.  We  do 
not  derive  the  difference  equations  here  because  they  can 
be  found  in  many  papers,  including  those  given  above. 
Instead,  we  consider  some  of  the  attractive  features  of 
FDTD  that  are  useful  for  our  hybrid  method. 

FDTD  models  the  actual  real-time  behavior  of  EM 
fields  which  makes  it  suitable  for  impulsive  and  transient 
analysis.  The  frequency  spectrum  can  also  be  easily  ob¬ 
tained  via  Fourier  transformation  of  the  FDTD  time 
solution.  In  contrast,  the  use  of  frequency-domain  meth¬ 
ods  to  get  the  frequency  spectrum  or  time  variation  re¬ 
quires  solving  the  problem  repeatedly  for  multiple  fre¬ 
quencies  in  the  frequency  range  of  interest. 

In  FDTD,  the  field  computations  involve  only  simple 
arithmetic  operators  with  no  matrix  manipulation.  The 
number  of  floating-point  operations  per  time  step  and  the 
memory  storage  are  both  proportional  to  the  number  of 
unknowns,  n,  compared  to  and  respectively,  in 
frequency-domain  methods  such  as  the  method  of  mo¬ 
ments  (assuming  LU  decomposition).  The  equations  in 
FDTD  can  also  be  easily  vectorized  for  the  Cray  com¬ 
puter.  Despite  these  advantages,  FDTD  is  computation¬ 
ally  costly  for  solving  electrically  large  problems.  For 
example,  a  rectangular  cavity,  with  a  cross  section  of 
30A  X  40A  and  a  depth  (or  length)  of  lOOA,  requires  5.76 
billion  unknowns  for  a  mesh  density  of  20  cells/A.  More¬ 
over,  the  FDTD  field  equations  have  to  be  repeatedly 
evaluated  for  at  least  8000  time  steps  (or  iterations,  as¬ 
suming  two  time  steps  for  the  field  to  traverse  from  one 
cell  to  another  along  the  length  of  the  cavity)  before  the 
interior  irradiated  field  can  be  observed  at  the  aperture  of 
the  cavity.  Finally,  FDTD  cannot  easily  model  curve 
boundaries/surfaces  although  some  solutions  have  been 
proposed  to  partly  overcome  this  shortcoming  [16],  [17]. 

C.  Coupling  between  Ray  and  FDTD 

As  was  alluded  to  earlier,  FDTD  can  be  computation¬ 
ally  inefficient  when  it  is  used  to  solve  electrically  large 
geometries.  For  the  previous  example  of  the  rectangular 
cavity,  most  of  the  computational  time  and  effort  are 
wasted  in  the  two-way  propagation  of  the  field  along  the 
length  of  the  cavity.  The  proposed  hybrid  ray-FDTD 
method  overcomes  this  inefficiency  by  using  a  ray  method, 
specifically  GRE,  to  propagate  the  fields  from  the  aper¬ 
ture  to  the  (vicinity  of  the)  termination.  However,  since  it 
is  difficult  to  apply  ray  methods  to  complex  terminations, 
FDTD  is  used  to  determine  the  interaction  of  the  cavity 
fields  with  the  termination. 

As  with  most  EM  hybrid  methods,  the  coupling  be¬ 
tween  the  individual  methods  has  to  be  handled  with  care. 
Consider  the  cavity  shown  in  Fig.  2.  The  cavity  has  been 
divided  into  three  regions  corresponding  roughly  to  the 
air  intake,  the  engine,  and  the  exhaust  sections  of  a  jet 
engine  inlet.  The  imaginary  surface  Sj^  separates  regions 
1  and  2,  while  5^^  separates  regions  2  and  3.  In  region  1, 
the  cavity  is  assumed  to  be  smoothly  varying  for  high- 
frequency  methods  like  ray  methods  to  be  valid.  Ray 


tubes  are  traced  from  the  aperture  of  the  cavity  to  Sj, 
where  they  are  summed  to  form  a  high-frequency  solution 
of  the  cavity  fields  across  Sj^,  This  solution  on  is  used 
as  the  excitation  for  the  FDTD  computation  in  region  2. 

We  will  now  describe  two  approaches  to  the  coupling 
between  the  GRE  and  the  EDIT)  methods.  Depending  on 
the  size  of  a  ray  tube  when  it  reaches  Sy^,  it  will  either  not 
intersect  any  of  the  FDTD  grid  points  on  Sy^  or  it  will 
intersect  one  or  more  of  those  grid  points.  Therefore, 
different  ray  tubes  will  contribute  differently  to  the  total 
field  of  a  grid  point.  Moreover,  it  is  possible  that  some 
grid  points  will  not  have  any  ray  tubes  intersecting  them. 
As  a  result,  the  incident  fields  on  Sy  evaluated  via  rays 
will  not  be  smoothly  varying.  In  order  to  obtain  a  smoothly 
varying  as  well  as  an  accurate  field  on  Sy^^  the  size  of  the 
ray  tube  can  be  restricted  so  that  it  will  only  intersect  one 
FDTD  grid  point.  Alternatively,  some  form  of  interpola¬ 
tion  can  be  applied  to  the  ray  tubes  that  intersect  more 
than  one  grid  point.  In  both  of  these  approaches,  a  ray 
tube  is  launched  and  tracked  via  its  central  ray  to  Sy .  The 
projected  ray  tube  area  (for  the  3D  case;  width  for  the  2D 
case)  on  Sy^  is  then  determined.  If  this  area  is  greater 
than  some  specified  area,  the  ray  tube  is  subdivided 
and  the  process  of  tracking  and  determining  the  projeaed 
area  of  the  smaller  ray  tubes  is  repeated.  In  the  first 
approach,  is  equal  to  the  FDTD  grid  spacing  (usually 
A/20)  for  the  2D  case  so  that  if  the  projected  ray  tube 
width  is  less  than  Ap,  that  ray  tube  can  intersect  at  most 
one  FDTD  grid  point  on  Sy^.  In  such  an  event,  the  field 
specified  by  the  central  ray  will  be  added  to  the  inter¬ 
sected  grid  point.  Ray  tubes  that  do  not  intersect  any  grid 
point  are  ignored. 

For  the  second  approach,  Ap  is  larger,  but  no  greater 
than  A/2  and  (A/2)^  for  the  2D  and  3D  cases,  respec¬ 
tively,  for  reasons  given  in  [2].  In  our  2D  implementation, 
we  have  used  Ap  =  A/4  for  greater  accuracy  in  the 
evaluation  of  the  fields  on  Sy^.  Since  a  larger  ray  tube  may 
intersect  more  than  one  FDTD  grid  point  on  Sy,  its 
contribution  to  the  fiel$l  on  Sy  is  determined  by  convert¬ 
ing  its  ray  field  into  modal  fields.  Specifically,  the  contri¬ 
bution  of  the  ray  field  to  the  modal  coefficients  are 
determined  by  integrating  the  ray  field  over  the  projected 
area  of  the  ray  tube  assuming  a  linear  phase  variation  in 
the  field  over  the  projected  area  with  respect  to  the  field 
of  the  central  ray.  This  approach  assumes  that  the  fields 
on  Sy^  are  expressible  in  terms  of  parallel-plate  wave¬ 
guide  modes  for  2D  problems.  This  assumption  can  usu- 
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ally  be  satisfied  by  a  suitable  choice  of  S^-^.  In  realistic  3D 
problems,  there  is  usually  a  narrow  section  in  front  of  the 
termination  which  is  cylindrical  so  that  the  fields  in  this 
narrow  section  can  also  be  expressed  in  terms  of  modes. 
In  any  case,  when  the  modal  coefficients  have  been  ob¬ 
tained  by  summing  the  contributions  due  to  all  the  ray 
tubes,  the  desired  field  at  each  FDTD  grid  point  on 
can  be  determined. 

Comparing  the  two  approaches,  it  is  clear  that  the  first 
approach  requires  more  ray  tracing  (which  means  more 
computational  time  and  storage)  than  the  second  ap¬ 
proach  because  of  the  smaller  Ap.  However,  once  the  rays 
are  traced,  the  first  approach  uses  only  a  simple  summa¬ 
tion  of  the  ray  fields  to  obtain  the  desired  fields  at  the 
FDTD  grid  points  on  Sj^.  The  second  approach  is  less 
efficient  in  this  latter  aspect  because  it  has  to  compute  the 
modal  coefficients  and  then  sum  the  modal  fields  to  get 
the  desired  fields.  However,  in  realistic  3D  problems,  the 
first  approach  may  not  be  viable  because  the  amount  of 
ray  tracing  can  become  overwhelmingly  excessive. 

Another  consideration  in  the  coupling  of  the  GRE  and 
FDTD  methods  is  the  selection  of  a  suitable  time  varia¬ 
tion  for  the  excitation  since  the  former  is  a  frequency- 
domain  method  while  the  latter  is  a  time-domain  method. 
There  are  two  possible  schemes  for  the  time  variation:  the 
sinusoidal  steady-state  time  variation  [18]  and  the  pulsed 
(usually  Gaussian  or  raised-cosine)  time  variation.  For  the 
steady-state  FDTD,  the  ray  solution  on  Sf^  is  evaluated 
only  at  a  single  frequency  of  interest.  Based  on  the 
complex  ray  field  solution  at  Sf^,  the  excitation  can  be 
made  to  vary  sinusoidally  with  time.  For  the  pulsed  FDTD, 
there  are  two  possible  alternatives.  For  the  first  alterna¬ 
tive,  the  ray  solution  on  Sr^  is  computed  over  a  range  of 
frequencies  corresponding  to  the  frequency  content  of  the 
pulse.  The  excitation  can  then  be  obtained  by  an  inverse 
Fourier  transform  of  the  product  of  the  ray  solution  and 
the  Fourier  transform  of  the  pulse.  This  alternative  is  not 
attractive  because  the  resultant  inverse  transform  will 
have  a  wide  time  window  with  a  number  of  significant 
pulses  due  to  the  different  arrival  times  at  of  the 
reflected,  diffracted,  and  reflected-diffraction  fields.  A 
better  alternative  is  to  use  a  basis  (e.g.,  modes)  as  excita¬ 
tion  for  the  pulsed  FDTD  to  characterize  the  termination 
section  in  terms  of  a  termination  scattering  matrix.  This 
scattering  matruc,  together  with  the  ray  solution  on  57- 
(expressed  also  in  terms  of  the  basis  set),  can  then  be  used 
to  find  the  cavity  scattered  field.  The  pulsed  time  variation 
scheme  is  more  efficient  for  problems  which  •require  mul¬ 
tiple  frequency  solutions  while  the  steady-state  scheme  is 
more  efficient  for  problems  which  require  only  a  single 
frequency  solution. 

Regardless  of  the  time  variation  used,  the  excitation 
produces  a  wave  which  propagates  toward  the  termination 
and  interacts  with  it.  For  the  geometry  shown  in  Fig.  2, 
part  of  the  wave  may  be  transmitted  to  region  3  through 
Sj^  while  the  remainder  is  reflected  back  toward  Sj^.  If  we 
assume  that  the  waves  leaving  region  2  through  the  imagi¬ 
nary  surfaces  and  Sy^  do  not  return,  then  an  absorb¬ 


ing  boundary  condition  (ABC)  such  as  the  ones  intro¬ 
duced  by  Higdon  [19],  [20]  or  Mur  [21]  can  be  applied  in 
the  FDTD  computations  at  each  of  the  two  surfaces.  The 
above  assumption  is  reasonable  because  most  jet  engine 
inlets  are  shaped  in  such  a  way  that  there  is  very  little 
energy  that  returns  to  region  2  upon  its  exit  from  there. 
Otherwise,  we"  can  convert  the  waves  leaving  region  2 
back  into  rays  (using  GRE)  and  track  those  rays  that 
return  to  region  2.  These  returning  rays  act  as  an  addi¬ 
tional  excitation. 

D.  Absorbing  Boundary  Conditions  at  Sy^  and  Sy^ 

As  mentioned  previously,  absorbing  boundary  condi¬ 
tions  (ABCs)  are  applied  at  the  imaginary  boundary  sur¬ 
faces  Sy  and  Sy^  for  the  proper  transmission  of  waves 
through  these  surfaces.  However,  the  two  ABCs  cannot  be 
implemented  in  the  same  way  because  of  different  field 
conditions  at  their  respective  boundary  surfaces.  With 
reference  to  Fig.  2,  the  ABC  at  Sy^  has  to  properly 
account  for  the  waves  which  pass  through  Sy^  in  both 
directions;  it  has  to  account  for  the  incident  excitation  af 
Sy  and  the  scattered  field  due  to  the  termination.  In 
contrast,  the  absorbing  boundary  at  Sy^  only  has  waves 
transmitted  through  it  from  region  2-  to  region  3  assuming 
that  there  are  no  waves  transmitted  through  Sy,  from 
region  3  to  region  2  (recall  discussion  in  previous  section). 

The  ABC  at  has  to  transmit  (or  absorb)  the  waves 
that  are  scattered  by  the  termination  toward  Sy^  without 
destroying  or  affecting  the  incident  excitation  at  Sy^.  We 
will  demonstrate  how  this  function  can  be  accomplished 
with  the  second-order  ABC  of  Mur  [21]  in  a  2D  problem. 
Applying  the  ABC  given  by  (17)  of  [21]  to  the  scattered 
electric  field  ,  (for  the  TM^  case)  at  Sy^,  we  have 

£;.:'(0,/)  =  £;.,(!, ;)fl, [£;.:'(!, y)  -  £;.,(o,y)] 

-fl,[//;.:'/2(o,y  0  -//;:'/2(o,y  -  {) 

+  +  ^)  -  -  D] 

(4) 

(following  Mur’s  FDTD  indexing;  see  Fig.  3  also  for 
indexing),  where  E, ,  and  //,  ,  are  the  fields  scattered  by 
(or  reflected  from)  the  termination.  The  constants  a,  and 
a,  are  given  by 

Cq  At  -  Ax 

Co  At  -t-  Ax  ’ 

_  /JLpC^AtAx 

2  Ay  (cp  At  +  Ax) ' 

Replacing  the  scattered  field  components  (£j.,,  //,.,)  with 


(5) 
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Fig.  3.  Geometiy  of  a  wedge  with  respect  to  the  FDTD  grid  for  the  TE, 

case. 


where  the  subscript  /  in  denotes  the  incident 

component  of  the  total  field,  (4)  becomes 

•[£r'(l,;)-£;;'(l,;)-£;.,(0,;)] 

;  +  i) 

+  {}-  - 1) 

-  i)  +  -  \)]. 

(6) 

Except  for  £^'(1,;)  and  all  the  other  field 

components  in  (6)  are  computed  values  for  the  previous 
(one  or  half)  time  step.  £"‘^'(l,y)  can  be  calculated  for 
the  current  time  step  from  the  regular  (total-field)  differ¬ 
ence  equation.  £"|'(l,y)  can  be  obtained  via  the  ray 
method  (as  was  done  for  the  excitation  at  the  absorbing 
boundary)  or  from  the  propagation  of  the  incident  excita¬ 
tion.  After  £^ ,  on  has  been  found,  the  total  field 
components  of  and  in  region  2  can  be  updated  via 
the  regular  FDTD  equations. 

The  procedure  described  above  for  obtaining  the  scat¬ 
tered  fields  at  due  to  the  reflection  of  the  fields  from 
the  termination  has  to  be  used  whenever  the  steady-state 
sinusoidal  time  variation  is  chosen  for  the  FDTD  compu¬ 
tation  in  region  2.  However,  this  procedure  may  or  may 
not  be  necessary  when  the  pulsed  time  variation  is  chosen, 
depending  on  the  time  window  of  the  excitation  and  the 
closest  distance  between  Sj^  and  the  termination. 


The  ABC  at  Sj^  does  not  require  any  spiecial  treatment 
like  the  one  at  if  we  assume  that  waves  are  transmit¬ 
ted  through  only  from  region  2  to  region  3  and  not 
vice  versa.  In  this  case,  the  transmitted  field  is  also  the 
total  field.  Therefore,  any  suitable  ABC  for  the  absorption 
of  the  total  field  at  can  be  applied.  Note  that  the  ABC 
(4)  can  also  be  applied  at  Sj  with  the  scattered  field 
variables  replaced  appropriately  by  the  total  field  vari¬ 
ables  (with  the  proper  spatial  indices). 

£.  Scattered  Field  Computation 

To  determine  the  cavity  scattered  fields,  the  appropri¬ 
ate  field  solutions  in  regions  1  and  2  have  to  be  used.  One 
possible  way  of  finding  the  cavity  scattered  field  is  to 
launch  GRE  rays  into  region  1  using  the  FDTD  solution 
at  Sj-^.  These  rays  are  tracked  to  the  front  aperture  of  the 
cavity  where  aperture  integration  can  be  applied  to  find 
the  scattered  field.  Unfortunately,  this  method  requires 
that  rays  be  traced  both  into  the  cavity  for  the  incident 
excitation  at  St^  and  out  of  the  cavity  for  the  cavity 
scattered  field  computation. 

A  more  suitable  and  efficient  way  for  finding  the  cavity 
scattered  field  is  based  on  the  termination  reciprocity 
integral  developed  by  Pathak  and  Burkholder  [22].  This 
integral  is  given  by 

E,(P)  •/»,£/  iE,XH,-E,X  H,)  •  h  ds,  (7) 

where  E,iP)  is  the  desired  cavity  scattered  field  at  the 
observation  point  P  and  P,  is  the  strength  of  an  electric 
current  point  (test)  source.  (E„H,)  are  the  fields  scat¬ 
tered  by  the  termination  in  the  cavity  while  (E,,H,)  are 
the  fields  radiated  by  the  test  source  in  the  presence  of 
the  cavity  structure  without  the  termination,  n  is  the  unit 
vector  normal  to  the  surface  as  shown  in  Fig.  2.  The 
approximation  in  (7)  assumes  that  the  source  and  ob¬ 
server  are  in  direct  view  of  the  open  front  end  so  that  the 
contribution  to  £,(£)  from  the  scattered  field  exiting 
through  is  negligible  compared  to  the  fields  exiting 
through  S-p. 

To  see  tlie  usefulness  of  (7),  consider  how  the  cavity 
scattered  fields  are  obtained  using  a  purely  ray-based 
approach.  For  example,  in  the  original  implementation  of 
SBR  [1],  [7],  two-way  (in  and  out  of  the  cavity)  ray  tracing 
has  to  be  performed  before  the  scattered  field  can  be 
determined  via  aperture  integration.  With  the  reciprocity 
integral,  the  rays  are  only  traced  from  the  open  end  of  the 
cavity  to  the  termination  and  back  to  Sp^.  Therefore,  the 
amount  of  ray  tracing  is  basically  reduced  by  about  half. 
This  reduction  is  even  more  significant  when  the  reciproc¬ 
ity  integral  is  used  in  the  hybrid  method  because  the  rays 
are  only  traced  from  the  open  end  to  Sp^. 

From  (7),  we  see  that  the  integral  is  independent  of  the 
method(s)  used  to  find  the  two  sets  of  tangential  fields  on 
Sp^.  It  is  particularly  easy  to  apply  the  integral  to  find  the 
backscattered  field  in  our  hybrid  method  since  both  sets 
of  fields  at  Sp^  are  readily  available.  Specifically,  for 
backscatter  computation,  iE„H,)  is  the  initial  excitation 


LEE  AND  CHIA:  HYBRID  RAY-FDTD  METHOD 


1565 


on  evaluated  by  the  GRE  method  for  the  FDTD 
computation  in  region  2;  its  evaluation  has  already  been 
described  earlier.  (£,,  ff,)  is  the  termination  scattered 
field  that  is  transmitted  through  from  region  2.  It  has 
to  be  evaluated  on  Sj-^  from  the  FDTD  algorithm.  How¬ 
ever,  the  tangential  E  and  H  fields  in  a  Yee  cell  lie  on 
different  planes.  In  particular,  the  tangential  E  and  H 
planes  are  half  a  cell  width  apart  in  the  direction  normal 
to  both  of  these  planes.  Consider  for  example,  the  2D  TE^ 
case  shown  in  Fig.  3,  where  the  tangential  magnetic  field 
H.(0,j)  lies  on  and  the  tangential  electric  field  E/j,;) 
lies  on  a  plane  parallel  to  at  Ax/2  away.  To  obtain 
,(0,y)  (which  is  symbolically,  the  termination  scattered 
field  j  on  57^),  an  extrapolation  of  the  neighboring 
values  of  ,  is  used.  Specifically,  we  assumed  that  the 
gradient  of  the  scattered  field  Ey  ,  in  the  normal  direction 
of  Sy^,  dE,  ,/dx,  at  X  =  Ax/4  is  equal  to  that  at  x  =  Ax, 
so  that 

EyJOJ)  =  K/i.y)  -  \EyJlj).  (8) 

For  the  steady-state  FDTD,  the  scattered  field  compo¬ 
nents  (£y  j)  on  the  right-hand  side  of  (8)  are  replaced  by 
the  difference  between  their  respective  total  (£^)  and 
incident  (£,,.,)  field  components.  This  extrapolation 
scheme  can  be  similarly  applied  to  the  tangential  electric 
or  magnetic  field  in  the  3D  case.  When  all  the  tangential 
fields  within  the  integral  of  (7)  have  been  appropriately 
manipulated  into  the  frequency  domain,  the  cavity  scat¬ 
tered  field  E,(P)  can  be  computed. 

In  summary,  the  hybrid  ray-FDTD  method  uses  high- 
frequency  ray  solutions  as  input  excitation  for  the  evalua¬ 
tion  of  the  termination  scattered  field  via  the  low- 
frequency  FDTD  method;  the  scattering  of  the  cavity  is 
then  obtained  via  the  termination  reciprocity  integral.  The 
hybrid  ray-FDTD  method  combines  the  efficiency  of  the 
ray  method  (and  the  reciprocity  integral)  with  the  model¬ 
ing  flexibility  of  the  FD1T)  method.  At  the  same  time,  it 
overcomes  some  of  the  limitations  of  the  individual  meth¬ 
ods  in  analyzing  the  scattering  from  a  cavity. 

III.  Results 

The  hybrid  GRE-FDTD  method  has  been  imple¬ 
mented  in  three  different  programs  (see  Table  I)  using 
various  combinations  of  ray  tracing  and  time  variations  to 
determine  the  radar  cross  section  (RCS)  of  parallel-plate 
cavities  with  different  PEC  plug  terminations.  Program  1 
implements  the  GRE  method  such  that  the  projected 
areas  of  the  ray  tubes  are  less  than  A/20  to  find  the 
incident  field  on  Sy^  while  the  other  two  programs  imple¬ 
ment  the  GRE  method  such  that  the  ray  fields  are  con¬ 
verted  into  modes  to  find  the  same  fields.  Programs  1  and 
2  implement  the  FDTD  algorithm  using  the  steady-state 
sinusoidal  time  variation,  while  program  3  implements  the 
FDTD  algorithm  using  the  basis-pulsed  time  variation 
approach. 

In  all  three  programs,  the  aperture  of  the  cavity  is 
divided  into  three  subapertures  from  which  rays  are 


TABLE  I 

Different  Implementations  of  the  Hybrid 
GRE-FDTD  Method 


Program 

GRE 

FDTD 

1 

A.  -  A/20 

sinusoidal 

2 

Ay  -  A/4 

sinusoidal 

3 

Ay  -  A/4 

pulsed-basis 

launched  within  an  angle  of  ±75°  of  the  cavity  axis  to 
compute  the  excitation  at  Sy^.  Sy^  is  fixed  at  4A  from  the 
closed  end  of  the  cavity.  For  the  FDTD  algorithm,  a 
A/20  X  A/20  spatial  grid  is  used  (unless  it  is  specified 
otherwise)  with  a  time  step  of  Ax/(2co);  Mur’s  [21]  ABC 
is  applied  at  Sy^  to  absorb  the  fields  scattered  toward  Sy^ 
by  the  termination.  No  ABC  is  implemented  at  Sy^  since 
region  3  (see  Fig.  2)  is  nonexistent  for  our  terminated 
cavity.  The  scattered  field  is  obtained  via  the  reciprocity 
integral  (7)  which  is  now  exact  since  is  perfectly 
conducting  for  the  terminated  (PEC)  parallel-plate  wave¬ 
guide. 

In  the  following  examples,  we  will  only  present  results 
obtained  via  program  2.  These  results  are  compared  with 
reference  solutions  obtained  via  the  hybrid  asymptotic 
modal-method  of  moments  (modal-MM)  [23].  The  results 
from  the  other  two  programs  are  not  presented  here 
because  there  is  little  difference  between  their  results  and 
that  of  program  2.  In  fact,  the  purpose  of  programs  1  and 
3  is  to  check  the  accuracy  of  the  approach  employed  in 
program  2.  Moreover,  program  3  validates  the  time- 
domain  (pulsed-basis)  approach  which  has  great  potential 
in  (2D  and  3D)  scattering  problems  where  the  cavity  RCS 
is  required  at  multiple  frequencies.  For  the  purpose  of 
this  paper,  the  results  from  program  2  are  sufficient  to 
illustrate  the  utility  and  accuracy  of  the  hybrid  method. 

Figures  4  and  5  show  the  RCS  patterns  of  a  9.6A-wide 
and  30A-long  parallel-plate  waveguide  cavity  with  a  2A  x 
4A  rectangular  PEC  plug  termination  for  the  TM^  and 
TE^  cases,  respectively.  The  GRE-FDTD  patterns  in 
both  figures  show  excellent  agreement  with  the  corre¬ 
sponding  modal-MM  patterns.  The  slight  differences  be¬ 
tween  the  GRE-FDTD  and  the  modal-MM  solutions  can 
be  attributed  to  two  causes.  One  of  the  causes  is  the 
inability  of  the  ABC  to  absorb  the  higher-order  modes 
excited  by  the  rectangular  plug  termination.  These 
higher-order  modes  result  from  the  reflection  and  edge 
diffraction  of  the  fields  from  the  termination.  They  are 
particularly  significant  at  larger  incident  angles  since  the 
termination  reflected  fields  (which  are  the  incident  fields 
after  undergoing  multiple  reflections  in  the  termination 
section)  are  incident  at  large  angles  from  the  normal  on- 
the  absorbing  boundary  Sy^.  These  higher-order  modes 
are  not  transmitted  (or  absorbed)  as  well  as  the  lower- 
order  modes  by  the  second-order  Mur’s  ABC  (see  Table 
II  [24]).  In  addition,  their  multiple  reflections  between  the 
absorbing  boundary  and  the  termination  result  in  appre¬ 
ciable  errors  in  the  scattered  field.  Another  source  of 
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subapertures). 


Fig.  5.  Radar  cross-section  patterns  of  parallel-plate  cavity  with  a 
reaangular-shaped  PEC  plug  termination  using  the  hybrid  modal-MM, 
GRE-FDTD,  and  GRE  methods,  a  -  9.6A,  /  »  30A,  <f  -  4A,  •  2A, 

w  -  4A;  TE,:  — ,  modal-MM;  -n-  +,  GRE-FDTD;  — .  GRE  (three 
subapertures). 


error  is  the  ray  solution  on  Sj-^.  These  solutions  are  less 
accnirate  at  large  incident  angle  6  than  at  small  6.  This 
difference  in  accuracy  is  a  result  of  limiting  the  ray 
solution  to  those  rays  whose  launch  angles  are  within 
±75°  of  the  cavity  axis.  At  large  6,  this  (artificial)  angular 
limit  may  not  be  adequate.  For  example,  for  an  incident 
angle  of  60°,  the  magnitude  of  the  radiation  power  pattern 
of  each  subaperture  at  75°  is  -8  dB  (40%)  with  respect  to 
the  beam  maximum,  while  the  magnitude  of  the  first  and 
second  side  lobes  are  -9  dB  (35%)  and  -13  dB  (22%), 
respectively.  Therefore,  significant  contribution  by  rays 
whose  launch  angles  are  greater  than  75°  are  omitted 
from  the  final  ray  solution  on  Sj- . 

Also  shown  in  Figs.  4  and  5  are  the  respective  patterns 
obtained  using  the  GRE  method  alone.  In  this  purely 


TABLE  II 

Refixctjon  Coefficients  for  the  Mur  Absorbing  Boundary 
Condition  at  the  Modal  Angles  of  the  Higher-Order 
Modes  in  a  9.6A-Wide  Parallel-Plate 
Waveguide  Ovtty 


Approx,  modal  angle  of 
hi^er-order  modes  (deg) 

Approx,  refleaion 
coefficient  (percent) 

47 

3 

51 

5 

56 

8 

"  62 

10 

70 

25 

82 

50 

GRE  approach,  the  rays  are  launched  from  the  aperture 
and  traced  to  S-j-^  where  they  form  the  incident  field.  The 
rays  are  then  traced  beyond  S-j-^  to  the  termination  where 
they  are  scattered  back  to  Sj-^  to  form  the  termination 
scattered  field.  Both  the  incident  and  scattered  fields  are 
converted  into  modal  fields  before  they  are  used  in  the 
reciprocity  integral  to  obtain  the  cavity  scattered  fields.  It 
is  clear  from  Figs.  4  and  5  that  the  GRE  solutions  agree 
very  well  with  the  reference  and  hybrid  solutions  for  6  up 
to  about  15°.  This  agreement  suggests  that  the  diffracted 
fields  from  the  termination  are  negligible  for  small  6. 
Beyond  that,  the  GRE  solutions  show  significant  devia¬ 
tions  from  the  reference  solutions  at  certain  angles.  From 
these  results,  it  is  envisaged  that  for  a  complex  termina¬ 
tion,  the  hybrid  GRE-FDTD  method  will  produce  far 
more  accurate  results  than  the  purely  GRE  approach 
which  does  not  account  for  the  termination  diffracted 
fields. 

Figure  6  shows  the  resultant  RCS  patterns  for  the  TM, 
case  when  the  rectangular  plug  termination  is  replaced  by 
a  wedge-shaped  plug  termination.  The  height  and  base 
width  of  the  wedge-shaped  plug  are  2  A  and  4A,  respec¬ 
tively.  The  wedge  is  modeled  in  the  FDTD  code  with  the 
electric  field  tangent  to  the  wedge  surface.  The 
GRE-FDTD  pattern  is  again  in  excellent  agreement  with 
the  modal-MM  pattern. 

Figure  7  shows  the  RCS  patterns  of  the  wedge- 
terminated  cavity  for  the  TEj  case.  The  wedge  is  now 
modeled  in  the  FDTD  code  with  the  magnetic  field 
tangent  to  the  wedge  surface.  In  contrast  to  the  TM^  case 
where  the  tangential  electric  field  vanishes  on  the 
wedge  surface,  the  tangential  magnetic  field  for  the 
TEj  case  does  not  vanish.  Instead,  it  has  to  be  specially 
treated.  Specifically,  we  apply  Stoke’s  theorem  to  one  of 
Maxwell’s  curl  equations  to  obtain  a  difference  equation 
for  updating  the  tangential  correctly.  Applying  Stoke’s 
theorem  to  Maxwell’s  curl-E  equation  results  in 

r  '  r  dH 

®  Edl^-fioj  ~-2ds,  (9) 

^ABCA  ■'44BC  ot 

where  the  line  and  surface  integrals  are  taken  around  the 
closed  loop  ABCA  and  the  triangle  ABC,  respectively,  as 
shown  in  Fig.  3.  After  integration  and  rearrangement,  (9) 


becomes 


+CF.  (10) 


where  the  correction  factor  CF  is  given  by 


Notice  that  without  CF,  (10)  is  just  the  regular  FDTD 
equation  for  updating  Therefore,  CT  adjusts  for  the 
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e (deg) 


Fig.  10.  Radar  cross-section  patterns  of  parallel-plate  cavity  withy  a 
semicircular-shaped  PEC  plug  termination  using  the  hybrid  modal-MM 
and  GRE-FDTD  methods,  a  -  9.6A,  /  =  30A,  d  «  4A,  r  -  2A;  TE  :  — , 
modal-MM;  •+  -H-,  GRE-FDTD  (A/20);  OOO,  GRE-FDTD  (A/40) 
(three  subapertures). 


H/s  that  are  tangent  to  the  wedge  surface.  Once  again, 
the  GRE-FDTD  pattern  in  Fig.  7  shows  excellent  agree¬ 
ment  with  the  modal-MM  pattern  except  for  large  6. 

For  the  final  example,  a  PEC  semicircular  plug  of 
radius  2  A  is  substituted  for  the  wedge-shaped  plug.  The 
curved  surface  of  the  semicircular  plug  is  modeled  using 
its  stair-stepped  approximation  as  shown  in  Fig.  8  for  the 
TMj  case.  A  similar  approximation  of  the  semicircular 
plug  is  used  for  the  TE^  case.  The  RCS  patterns  of  the 
new  cavity  are  shown  in  Figs.  9  and  10  for  the  TM^  and 
TEj  cases,  respectively.  In  these  figures,  we  have  also 
plotted  the  patterns  obtained  with  the  smaller  A/40  FDTD 
grid.  The  patterns  obtained  with  the  A/20  grid  show  some 
semblance  to  the  reference  patterns,  while  those  obtained 
with  the  finer  grid  show  better  agreement  with  the  refer¬ 
ence  patterns,  particularly  for  the  TE^  case.  The  signifi¬ 
cant  improvement  in  the  latter  case  is  attributed  to  the 
fact  that  in  2D  geometries,  TE^  modes  suffer  dispersion 
due  to  the  stair-stepped  approximation,  while  TM^  modes 
does  not  [25].  This  dispersion  (or  numerical  error)  reduces 
when  a  finer  FDTD  grid  is  used.  Therefore,  the  pattern 
for  the  TEj  case  shows  marked  improvement,  while  that 
of  the  TM^  case  is  basically  unchanged. 

IV.  Summary 

In  this  paper,  we  have  introduced  a  hybrid  method 
which  combines  ray  methods  with  FDTD.  The  use  of  the 
hybrid  method  allows  us  to  overcome  the  deficiencies  of  a 
single  method.  The  FDTD  calculation,  which  is  computa¬ 
tionally  intensive,  is  limited  to  a  small  region  around  the 
termination.  The  ray  solution  is  used  to  evaluate  the  fields 
in  the  remainder  of  the  cavity  where  it  is  assumed  that  the 
geometry  is  simple  enough  for  high-frequency  approxima¬ 
tions  to  be  valid.  Results  have  been  presented  and  com¬ 


pared  to  a  hybrid  modal-MM  solution  to  demonstrate  its 
accuracy.  From  these  results,  it  is  evident  that  the  hybrid 
ray-FDTD  method  shows  great  potential  for  eventually 
providing  an  accurate  solution  for  the  electromagnetic 
scattering  from  a  realistic  three-dimensional  jet  engine. 
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ON  THE  CHOICE  OF  METRONS  IN  THE  MEI  METHOD 
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Department  of  Electrical  Engineering 
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I.  Introduction 

Finite  element  and  finite  difference  methods  find  increasing  application  to  the  numerical 
analysis  of  problems  involving  unbounded  solution  domains.  This  problems  require  special 
techniques  for  boundary  truncation.  A  new  way  to  truncate  the  boundary  has  been 
presented  in  [1]  through  the  use  of  the  measured  equations  of  invariance  (MEI).  Some  of 
the  theoretical  aspects  and  modifications  of  the  MEI  method,  as  well  as  several  examples 
of  its  practical  implementation,  have  been  reported  in  the  articles  [2-6]. 

We  shall  present  our  findings  and  the  fundamentals  of  the  method  using  the  example  of 
the  two  dimensional  electromagnetic  wave  scattering  by  a  perfectly  conducting  cylinder. 
We  shall  assume  TM  polarization,  as  shown  in  Figure  1.  A  general  linear  boundary 
condition  can  be  expressed  in  the  form; 

x:o.£f=o  (1) 

1=) 

where  El  stands  for  the  scattered  field  at  node  i.  and  M  is  the  number  of  nodes  coupled 
through  this  equation. 

What  is  specific  to  MEI  method  is  how  the  coefficients  o,  are  chosen.  Namely,  one  uses 
the  fact  that; 

£!  =  /  G,{s)Ms)ds,  0  <  s  <  1  (2) 

where  s  is  the  normalized  coordinate  along  the  conductor  surface,  Js(s)  is  the  surface 
current  induced  on  the  cylinder  and  G’,(#)  is  the  Green’s  function  from  point  s  on  the 
cylinder  to  node  i.  One  now  chooses  coefficients  o,  such  that  (1)  is  exactly  satisfied  for 
several  different  surface  currents  called  metrons; 

M  . 

Ta.  G,{s)M,(s)ds  =  0,  j  =  2 . M  (3) 

JdSc 
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We  need  M  —  \  metrons  since  one  of  the  MEl  coefficients  is  arbitrary  due  to  the 
homogeneous  nature  of  (1). 


II.  Null  Field  and  Its  Minimization 

We  have  shown  in  [2]  that,  although  (1 )  is  exactly  satisfied  for  metrons  A/j(s),  j  =  2, M , 
the  left-hand  side  of  (1)  is  not  equal  to  zero  for  all  possible  currents  J,(s).  We  should 
rather  write: 


Y,q,ei  =  r  (4) 

i=l 

where  R  is  the  residual.  Indeed,  by  using  (2),  (1)  becomes: 

[  EnMs)J.U)ds  =  R  (5) 

JaSc 

where: 

M 

E„uu{s)  =  '£^'GM  (6) 

t=j 

We  have  introduced  the  quantity  Enuii  in  (2]-  We  call  it  “the  null  field”  since  it  physically 
represents  the  incident  field  on  the  conductor  surface  due  to  M  line  sources  of  strengths 
Q,  placed  at  corresponding  nodes  and  because  E„uii  should  preferably  be  zero  if  (1)  is  to 
be  true  for  all  possible  currents  J,(s). 

Namely,  if  (1)  is  to  be  an  accurate  approximation  of  (4),  the  MEI  coefficients  should  be 
chosen  so  that  R  in  (5)  is  as  small  as  possible  for  all  possible  currents.  One  solution  has 
been  proposed  in  [2]  where  we  have  minimized  the  null  filed  in  the  least  square  sense: 


ds  =  min 


Using  (6),  this  condition  yields: 


(7) 


M  . 

j  G.(s)Gj(s)*ds  =  0,  j  =  2,...,M  (8) 

1=1  ■'SSc 

A  comparison  with  (-3)  reveals  that  this  is  the  same  result  as  would  be  obtained  by  using 
the  conjugate  of  Green’s  functions  as  metrons: 


M,{s)  =  G,(s)*. 


i  =  2 . M 


(9) 
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III.  Spectral  Approach  to  Null  Field  Minimization 

Condition  (7)  minimizes  residual  R  in  (5)  for  all  possible  currents  J,(s).  This  would  be 
the  optimum  choice  if  we  knew  absolutely  nothing  about  the  actual  current  induced  on 
the  cylinder.  In  terms  of  the  spectral  representation: 

+  0O 

=  E  0  <  S  <  1  (10) 

ns=  — oo 

this  assumption  would  mean  that  (7)  is  the  optimum  choice  when  all  components  of  the 
Fourier  expansion  in  (10)  are  present  in  the  actual  current. 

However,  this  is  not  the  case  in  general,  since  spectral  components  of  J,  whose  spatial 
variations  are  much  faster  than  that  of  the  incident  wave  are  usually  negligible.  It  may 
therefore  be  possible  to  devise  a  less  stringent  condition  than  that  in  (7)  and  obtain 
coefficients  that  would  be  more  suitable  for  a  particular  problem. 

This  can  be  done  as  follows.  Introduce  Fourier  expansions  of  Green’s  functions  and  null 
field: 


G.{s) 


+  00 


nss— oc 


(H) 


and: 


£nu,/(.o)=  E  (12) 

nsr—oo 

Equivalents  of  the  equations  (5)  and  (6)  in  the  spectral  representation  are  then: 

E  EnuUn^Scn)  =  ^ 

n=-oc- 


and: 


EnuU„  =  Eo.Gi„  (14) 

•=1 

Assuming  that  an  estimate  of  the  magnitudes  of  the  spectral  components  of  current  is 
available,  minimization  of  the  residual  R  in  (13)  in  the  least  square  sense  would  be: 


which  is  the  required  substitute  for  (7). 
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If  we  actually  perform  this  minimization  through  the  use  of  (14),  we  obtain: 

M  +O0 

=  0,  (16) 

1=1  n=-oo 

and  with  use  of  (11) 

M  .  +00 

5:q.  /  G.(s)  Y.  =  0,  j  =  2 . M  (17) 

JdSc  n=-oo 

Comparison  with  (3)  revelas  that  this  is  the  same  result  as  if  the  metrons  defined  by: 

+  00 

A/,(s)  =  (  5:  GjJJs,.„jV^'-)-,  j  =  2,...,M  (18) 

n=— oo 

were  used  to  find  the  coefficients. 

It  is  interesting  to  note  that  if  lJs„|  =  1  for  all  n.  indicating  the  presence  of  all  spectral 
components,  (18)  reduces  to  (9)  as  might  be  expected. 

On  the  other  hand,  if  IJgnI  =  1  for  |ij|  <  A'  and  |Jgn|  =  0  otherwise,  indicating  that 
high  spatial-frequency  components  are  negligible,  (18)  gives  exactly  the  same  results  as  if 
sinusoidal  metrons  up  to  A'"'  order  were  used: 

A/n(5)  =  n  =  -A . A  (19) 

and  the  overspecified  system  of  equations  for  coefficients  (3)  was  solved  using  the  least 
square  method. 


rv.  Numerical  Results 

We  have  encountered  three  different  types  of  metrons  in  the  preceeding  analysis.  Sinu¬ 
soidal  metrons  (19),  metrons  (9),  and  optimized  metrons  (18).  When  implementing 
optimized  metrons  it  is  more  computationaly  efficient  to  use  (16)  instead  of  (18). 

The  results  that  are  presented  in  this  section  are  obtained  using  a  finite  difference  mesh 
that  has  been  truncated  two  cells  from  the  conductor  surface.  Six  nodes  are  coupled 
through  (1),  and  therefore  A/  =  6  as  shown  in  Figure  1.  Nodal  density  W2is  chosen  to  be 
20  nodes  per  wavelength  and  all  integrations  were  performed  using  ten  integration  points 
per  cell. 

Results  in  Figures  2  and  3  were  obtained  for  the  case  of  a  circular  cylinder  whose  diameter 
is  10  wavelengths.  Figure  2  is  a  plot  of  llio  absolute  error  in  the  solution  for  surface 
current  density  J,  when  different  metrons  are  used.  Sinusoidal  metrons  for  N=2  give 
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slightly  better  solution  than  g'  metrons.  The  best  results  are  obtained  with  the  use  of 
optimized  metrons. 

When  implementing  (16),  we  have  estimated  the  magnitude  of  the  current  spectrum  |J*„| 
by  using  the  physical  optics  solution.  The  actual  spectrum  of  the  current  is  shown  in 
Figure  3.  This  figure  also  shows  the  spectrum  of  the  null  field  for  all  three  types  of 
metrons.  The  spectrum  of  the  null  field  corresponding  to  g’  metrons  is  very  large  in 
the  frequency  domain  of  interest,  although  it  is  minimal  when  considered  over  the  whole 
frequency  range.  On  the  other  hand,  the  null  field  spectrum  of  the  optimized  metrons  is 
minimal  exactly  over  the  frequency  range  that  is  most  important. 

\\e  might  expect  that  the  g'  metrons  will  give  better  results  than  the  same  number  of 
sinusoidal  metrons  when  the  current  contains  strong  high  frequency  components.  Figure 
4  shows  the  error  in  the  solution  for  a  rectangular  cylinder  with  a  high  aspect  ratio  (5A 
by  0.1  A).  Forty  sinusoidal  metrons  were  necessary  to  achieve  the  same  accuracy  as  with 
the  Sp"  metrons. 
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Figure  1.  Geometry  for  a  PEC  cylinder 
illuminated  with  a  TM  polarized  plane 
wave. 


sinusoidal  .V  =  2 


Figure  2.  Plot  of  the  error  in  the  sur¬ 
face  current  on  circular  cylinder  of  di¬ 
ameter  lOA.  Sinusoidal,  g’  and  opti¬ 
mized  metrons  are  compared.  Wave  is 
incident  normaly  at  the  point  s  —  0.5. 
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sinusoidal  .V  =  2 
g’  metrons 
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Figure  3.  Plot  of  the  spectra  of  null 
fields  corresponding  to  three  different 
sets  of  metrons  for  the  case  shown  in 
Figure  2.  The  normalized  spectrum  of 
the  surface  current  is  shown  for  refer- 


Figure  4.  Plot  of  the  error  in  the  sur¬ 
face  current  on  the  rectangular  cylinder 
5A  by  0.1  A.  Wave  is  incident  normally 
to  the  shorter  side.  Sinusoidal  metrons 
with  N  =  2  and  20  are  compared  with 
the  g'  metrons. 


ence. 


IEEE  TRANSACTIONS  ON  ANTENNAS  AND  PROPAGATION.  VOL.  42,  NO.  I.  AUGUST  1994 


1097 


A  Theoretical  and  Numerical  Analysis  of  the 
Measured  Equation  of  Invariance 

Jovan  O.  Jevti^,  Student  Member,  IEEE,  and  Robert  Lee,  Member,  IEEE 


Abstract — The  measured  equation  of  invariance  (MED  has 
been  recent!)  introduced  to  efhcientl)  and  accurately  handle  the 
boundary  truncation  for  finite  methods.  In  this  paper,  we  present 
a  theoretical  analysis  that  provides  several  important  insights  into 
the  capabilities  of  the  MEI.  From  the  numerical  study,  we  can 
explain  why  the  MEI  works  better  than  one  would  expect.  Both 
the  theoretical  and  the  numerical  analysis  demonstrate  that  the 
accuracy  of  the  solution  is  dependent  on  the  electrical  size  of 
the  geometry  as  well  as  the  distance  between  the  mesh  boundary 
and  the  geometry.  From  the  analysis,  we  propose  a  new  set  of 
metrons  that  is  less  sensitive  to  the  excitation  than  the  previously 
proposed  sinusoidal  metrons. 


I.  Introduction 

HE  finite  element  method  (FEM)  has  always  been  handi¬ 
capped  in  the  modeling  of  electromagnetic  phenomena 
because  of  its  inability  to  both  accurately  and  efficiently 
simulate  the  radiation  condition.  In  the  past  25  years,  re¬ 
searchers  have  developed  numerous  boundary  conditions  to 
model  the  radiation  condition.  The  work  can  be  separated 
into  two  classes  of  methods,  local  and  global.  In  the  local 
methods  [1]-15].  a  finite  difference  (FDj  equation  is  derived 
for  the  boundary  based  on  an  approximate  assumption  of 
the  behavior  of  the  field.  Usually,  the  wave  is  assumed  to 
be  outgoing  from  the  computation  domain  in  a  direction 
perpendicular  to  the  boundary  .  In  addition,  the  FD  equations 
are  derived  without  any  consideration  of  the  geometry  of 
interest.  The  major  advantages  of  the  local  methods  are  that 
they  maintain  the  sparsity  of  the  FEM  matrix  and  that  they  are 
easily  implemented  into  an  FEM  code.  The  major  disadvantage 
is  that  they  may  give  inaccurate  results  when  the  mesh 
boundary  is  placed  close  to  the  geometry.  To  assure  an  accurate 
solution,  the  boundary  must  be  extended  a  large  distance 
from  the  geometry  ,  which  may  dramatically  increase  both  the 
computation  time  and  the  storage  requirements.  The  choice  of 
where  to  place  the  boundary  is  also  not  clear.  Currently,  the 
only  way  to  determine  the  correct  placement  of  the  boundary 
for  complex  geometries  is  empirical.  The  minimum  distance 
between  the  scatterer  and  the  grid  boundary  is  dependent  upon 
the  geometry,  the  characteristics  of  the  excitation,  and  the 
degree  of  accuracy  desired;  so.  in  principle,  numerical  tests 
must  be  run  for  each  new  geometry. 
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The  global  methods  arc  based  on  a  rigorous  formulation 
in  which  all  the  nodes  on  the  boundary  are  coupled  together; 
therefore,  an  accurate  solution  can  be  obtained  without  regard 
to  the  distance  between  the  geometry  and  the  mesh  boundary  . 
In  some  methods  [6],  [7],  the  equations  that  couple  the  bound¬ 
ary  nodes  are  placed  directly  in  the  FEM  matrix.  However,  this 
has  the  effect  of  reducing  the  sparsity  of  the  matrix,  which  can 
result  in  a  significant  increase  in  computation  costs  In  other 
global  methods  [8],  [9],  the  boundary  equations  are  decoupled 
from  the  rest  of  the  FEM  matrix.  Although  these  methods 
are  usually  more  efficient  than  those  in  [6]  and  [7].  they  still 
require  the  solution  of  a  dense  matrix  due  to  the  fully  coupled 
boundary  equations.  In  addition,  the  computation  time  required 
to  decouple  the  boundary  equations  from  the  FEM  matrix  is 
comparable  to  the  solution  time  of  the  FEM  matrix  without 
the  boundary  equations. 

Recently,  the  measured  equation  of  invariance  (MEI)  [10] 
has  been  proposed  as  a  rigorous  local  method  for  modeling 
perfectly  conducting  cylinders.  It  is  based  on  the  concept  that 
accurate  local  FD  equations  can  be  written  at  the  boundary 
nodes  if  the  geometry  information  can  somehow  be  incor¬ 
porated  into  the  equations.  Thus,  a  general  FD  equation 
is  used  at  the  boundary  where  the  coefficients  in  the  FD 
equations  are  unknown.  These  coefficients  are  then  determined 
by  the  use  of  a  sequence  of  fields  called  measures.  A  more 
detailed  discussion  of  this  procedure  is  described  in  the  next 
section.  The  MEI  has  also  been  extended  to  dielectric  cy  linders 
Ill],  Although  both  groups  were  able  to  show  that  the  MEI 
works  well  for  several  geometries,  the  method  is  still  in 
its  infancy,  and  there  is  very  little  theoretical  foundation  to 
explain  why  the  MEI  works.  In  this  paper,  we  provide  some 
theoretical  insights  into  the  MEI.  We  show  that  the  MEI 
approach  still  produces  an  approximate  boundary  condition. 
However,  because  the  geometry  information  is  incorporated 
into  the  boundary  condition,  it  is  capable  of  producing  very 
good  results  with  appropriately  chosen  metrons.  The  sinusoidal 
meuons  turn  out  to  be  very  good  choices  for  the  geometry  we 
consider,  and  in  our  numerical  study  we  w  ill  explain  the  reason 
for  this.  The  numerical  study  also  demonstrates  the  effect  of 
changing  various  parameters  in  the  MEI. 

II.  Theoretical  Analysis 
A.  Revieys'  of  the  MEI 

To  analyze  the  MEI.  let  us  consider  the  problem  of  elec¬ 
tromagnetic  scattering  from  a  two-dimensional  perfectly  con¬ 
ducting  cylinder  in  free  space  (see  Fig.  1 ).  The  boundary  o1 
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the  cylinder  is  defined  to  be  dSr-  The  formulation  is  derived 
onl\  for  the  TM  (Ez.Hj.Hy)  polarization  since  the  same 
procedure  can  be  used  for  the  TE  case.  The  field  solution 
is  obtained  from  the  solution  of  the  Helmholtz  equation  for 
E-.  The  grid  surrounding  the  cylinder  can  either  be  a  finite- 
difference  or  a  finite-element  grid.  For  this  analysis,  let  us 
consider  a  finite-difference  grid.  At  the  boundary  of  the  grid, 
ue  cannot  apply  the  traditional  finite-difference  approximation 
because  the  traditional  difference  approximation  requires  grid 
points  that  are  exterior  to  the  grid.  Thus  we  must  find  another 
finite-difference  approximation  that  incorporates  the  radiation 
condition.  The  measured  equation  of  invariance  is  based  on 
the  principle  that  a  linear  equation  that  couples  the  boundary 
nodes  to  the  adjacent  nodes  can  be  formulated  without  using 
the  traditional  finite-difference  approximation.  One  possible 
coupling  is  the  six-node  coupling  shown  in  Fig.  1,  where  node 
6  is  coupled  to  nodes  I  through  5.  To  simplify  the  notation,  the 
analy  sis  is  presented  in  terms  of  this  coupling  scheme  without 
any  loss  of  generality.  The  linear  equation  at  node  6  can  be 
written  as 


=  o  (1) 

1=1 


where  £f,  is  the  scattered  electric  field  at  node  i  and  a,  are 
constant  coefficients  that  must  be  determined.  To  find  these 
coefficients.  Mei  et  al.  [10]  use  three  postulates  in  which  they 
conjecture  that  the  equation  in  ( 1 )  is 

1 )  location  dependent, 

2)  geometry  specific,  and 

3)  invariant  to  field  excitation. 

The  first  two  postulates  can  be  shown  to  be  true  by  numerical 
experimentation.  These  two  postulates  are  also  inherently 
incorporated  into  the  global  methods  described  in  the  previous 
section.  The  validity  of  the  third  postulate  is  not  as  clear. 
If  the  third  postulate  is  true,  then  the  coefficients  in  (1)  can 
be  determined  from  the  field  solutions  at  the  six  nodes  due 
to  five  linearly  independent  excitations  or  metrons  (only  five 
are  needed  since  the  sixth  coefficient  is  arbitrarily  specified), 
which  are  denoted  by  the  symbol  4'*.  The  five  metrons 
represent  five  linearly  independent  surface  current  density 
distributions  on  the  conducting  cylinder.  The  metrons  are 
usually  chosen  to  be  entire  domain  and  to  vary  sinusoidally 
over  the  cylinder.  For  example,  the  choice  in  [10]  is 

fcosTT/rs  fc  =  0,2,4  2 

^.sin  7r(A-  -f  l)s  =  1. 3 

where  0  <  s  <  1  describes  the  perimeter  of  the  cylinder. 
Although  the  results  presented  in  [10]  and  [11]  seem  to 
indicate  that  the  third  postulate  is  true,  there  were  some 
unexplainable  discrepancies  in  the  numerical  results.  The 
reasons  for  these  discrepancies  can  be  determined  from  an 
analysis  of  the  MEI. 


Fig.  I.  Perfectly  conducting  cylinder  geometry  for  analysis  of  the  MEI. 


B.  Analysis  of  the  MEI 

An  equation  can  be  written  relating  £*,  to  the  induced 
electric  surface  current  density  7,  on  the  cylinder,  as  follows; 

£*,  = -jWAto  /  9[Px-p')Jz{p')dt'  (3) 

JdS, 

where  p,  is  the  position  vector  to  node  i  and  9{p,.p')  is  the 
free-space  Green's  function,  which  is  given  by 

=  (4) 

The  variable  ko  is  the  free-space  wave  number  given  by  ko  ~ 
u)^fio(o.  The  above  equations  are  for  the  TM  polarization. 
For  the  TE  polarization,  the  nodal  unknown  is  the  z-directed 
magnetic  field,  which  is  related  to  the  electric  surface  current 
by  h  Vg  where  n  is  the  unit  normal  from  the  conductor 
surface.  Equation  (3)  is  substituted  into  (1)  to  obtain 

6  f 

53  a.  /  5(p..p')-/x(p')df'  =  0.  (5) 

t=i 


Recognizing  the  fact  that  g{p,,p')  =  g{p'.p,).  we  rewrite  (5) 
as 


a>9ip'’Pt) 


Mp')d('  =  0. 


(6) 


The  term  inside  the  brackets  has  a  special  physical  signifi¬ 
cance.  It  represents  the  electric  field  on  dS^  due  to  six  line 
sources  of  weights  a,  located  at  the  corresponding  nodes  i. 
Let  us  define  the  null  field  denoted  by  Enuii  to  be 
6 

^nuii(p')  =  p'€dSc.  (7) 

i=l 


In  order  for  the  third  postulate  to  be  true,  a  set  of  coefficients  a, 
must  be  found  such  that  (6)  is  true  for  all  possible  excitations. 
i.e.,  all  possible  Jz{p ')  for  p'  €  dSc-  This  is  possible  only 
in  the  instance  that 


Erxuwip')  =  0  p'  €  55,.  (8) 
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However,  it  is  impossible  to  choose  a  set  of  coefficients  o< 
such  that  (8)  is  satisfied.  Thus,  the  third  postulate  is  incorrect. 
]i  should  be  noted  that  we  use  the  notation  £nuii  in  (7)  since 
the  desired  value  of  Enuii  is  zero. 

Although  the  third  postulate  is  wrong,  the  numerical  results 
presented  in  [10]  indicate  that  the  method  has  some  merit.  The 
procedure  in  which  the  coefficients  are  determined  from  the 
metrons  may  provide  some  insight  into  the^accuracy  of  the 
MEl.  To  find  the  coefficients,  we  use  the  metrons  in  (2)  in 
place  of  J.  in  (5).  The  coefficient  are  then  determined  from 
the  resulting  matrix  equation.  With  this  choice  of  coefficients, 
the  null  field  is  orthogonal  to  any  current  distribution  on 
the  cylinder  that  can  be  represented  by  a  linear  combination 
of  the  metrons.  Thus,  any  excitation  of  the  cylinder  that 
induces  a  current  distribution  that  can  be  represented  by  a 
linear  combination  of  the  metrons  is  perfectly  absorbed  by  the 
boundary  condition  produced  from  the  MEl.  For  an  arbitrary 
case,  the  current  density  Jj  can  be  broken  down  into 

Jz  =  (9) 

where  is  the  part  of  the  current  density  that  can  be 
represented  by  the  metrons,  i.e., 

5 

Jme.  (10) 

fc=l 

Jjfi.  is  the  remainder  of  the  current  that  produces  the  error  in 
the  solution.  Once  the  coefficients  a,  have  been  determined, 
the  error  can  be  directly  incorporated  into  the  boundary 
equation  in  (1): 

6 

=  /?  (11) 

1=1 

where  the  residual  R  is  defined  from  (6)  and  (7)  to  be 

R  =  [  02) 

JdS, 

Since  is  orthogonal  to  £nuii-  02)  reduces  to 

R=  f  Ez,uu{p')J,.sd('.  03) 

Ji)S, 

For  small  smooth  cylinders.  J,e5  is  small  compared  to 
because  the  current  distribution  can  be  accurately  described  by 
the  metrons  in  (2).  Even  for  cases  where  is  locally  large, 
such  as  near  sharp  comers  of  conductors.  R  may  still  be  a 
small  value  since  the  integration  in  (13)  is  performed  over  the 
entire  cylinder.  For  cylinders  that  are  electrically  large.  J,es 
may  become  larger  than  .  which  may  result  in  significant 
errors  in  the  solution  unless  the  residual  R  due  to  J,r^  is  small. 

C.  Alternative  Metrons 

The  meuons  in  (2)  seem  to  have  been  chosen  arbitrarily 
because  of  the  assumptions  of  the  third  postulate.  This  was  a 
fortunate  choice,  as  we  will  see  in  the  next  section.  However, 
with  the  above  analysis  of  the  MEl.  it  is  now  possible  to 
develop  a  variety  of  different  schemes  for  the  determination 
of  the  coefficients  Ui.  In  this  section,  we  propose  a  scheme 


to  determine  the  coefficients  without  the  use  of  metrons.  The 
ideal  boundary  condition  would  be  the  case  where  coefficients 
can  be  found  such  that  (8)  is  satisfied.  Since  this  case  is  not 
possible,  one  choice  is  to  find  coefficients  such  that  fnun  is 
minimized  in  some  manner.  By  minimizing  £nuii  on  dSc,  the 
residual  £  in  (12)  is  minimized  independent  of  J^.  A  good 
method  for  the  minimization  of  £nuii  the  method  of  least 
squares,  in  which  the  coefficients  are  chosen  such  that 


df  =  Miniiiiuin. 


(14) 


Substituting  (7)  into  (14)  and  applying  the  method  of  least 
squares,  we  obtain  the  following  set  of  equations  for  o,; 


a,g(p'.p,)9’{p'.Pj)dC'  =  0 


j  =  I.2.-  -.o 
(15) 


where  g'  is  the  conjugate  of  the  free-space  Green's  function 
and  06  is  assumed  to  be  specified.  It  is  interesting  to  note  that 
the  determination  of  the  coefficients  from  (15)  is  mathemati¬ 
cally  equivalent  to  choosing  p*  as  the  metron.  Thus,  p*  can  be 
used  as  an  alternative  set  of  metrons  although  the  derivation  is 
based  on  the  minimization  of  £nuii  rather  than  on  any  concept 
of  “measuring"  the  FD  equations. 


Ill.  Numerical  STL'oy 

The  major  computation  cost  associated  with  the  MEl  is  the 
integration  of  the  metrons.  The  complexity  of  the  integral 
evaluation  is  SM  where  A'  is  the  number  of  nodes  on  the 
boundary  of  the  grid  and  M  is  the  number  of  integration  points 
on  the  cylinder  surface.  We  can  choose  M  —  A'  to  perform 
the  integration.  However,  the  coefficients  a,  are  very  sensitixe 
to  the  integral  evaluation,  and  we  found  that  with  the  choice  of 
JV/  =  A’  (20  nodes  per  wavelength),  the  results  are  inaccurate 
when  sinusoidal  metrons  are  used  and  the  number  of  layers 
is  less  than  four.  We  were  able  to  obtain  good  results  all  the 
way  down  to  two  layers  with  =  4.V.  The  inaccuracy  is  due 
to  the  singularity  of  the  Green's  function,  and  we  expect  that 
good  results  can  be  obtained  with  M  ft:  A'  if  unequally  spaced 
integration  points  are  used  or  if  the  singularity  is  extracted 
analytically  .  For  the  numerical  study  here,  our  interests  are 
not  in  numerical  efficiency  but  rather  numerical  accuracy; 
therefore,  we  have  chosen  M  =  lO.V  to  generate  the  resultN. 
The  MEl  is  not  a  symmetric  formulation,  so  the  resulting  finite 
difference  matrix  equation  is  also  nonsymmetrieal 

A.  Initial  Observations 

Numerical  problems  that  will  be  considered  in  this  paper 
all  relate  to  plane  wave  scattering  from  a  perfectly  conduct¬ 
ing  circular  cylinder  (see  Fig.  2).  This  geometry  is  chosen 
because  the  series  solution  is  known  and  also  because  the 
sinusoidal  metrons  are  easily  identified  with  terms  in  the  series 
solution.  We  will  always  use  a  nodal  density  of  20  nodes 
per  wavelength,  two  layers  of  nodes,  and  sinusoidal  metrons 
given  by  (2).  unless  explicitly  stated  otherwise.  It  should  be 
noted  that  the  conclusions  drawn  for  the  numerical  results  are 
not  applicable  to  problems  in  which  the  truncation  boundary 


Fif  Perfecil)  conducling  circular  cylinder, 

is  concave.  Up  to  this  point,  we  have  been  unsuccessful  in 
modeling  deep  cavities  with  the  MEI. 

The  application  of  the  MEI  method  to  small  smooth  cylin¬ 
ders  should  yield  excellent  results.  This  is  demonstrated  in  Fig. 

where  the  magnitude  of  surface  current  density,  computed 
u^ing  the  MEI  method,  is  compared  to  the  series  solution  for 
the  case  of  a  1-wavelength  diameter  circular  cylinder  where  the 
incident  plane  wave  is  chosen  to  be  TE  polarized.  Previously, 
we  argued  that  the  reason  for  the  good  agreement  lies  in 
the  fact  that  the  actual  current  distribution  can  accurately  be 
represented  h\  a  linear  combination  of  metrons;  therefore, 
there  should  be  very  little  error  since  the  MEI  boundary 
equations  work  perfectly  for  the  metrons.  To  verify  this 
fact  we  require  a  quantitative  measure  of  how  well  the 
MEI  boundary  equation  absorbs  field  radiated  by  a  particular 
current  distribution.  This  measure  is  provided  by  the  residual 
R.  which  has  been  introduced  in  (12).  A  smaller  residual 
indicates  a  more  accurate  boundary  condition  for  that  current 
distribution  on  the  cylinder.  Let  us  define  to  be  the 
harmonic  residual  due  to  an  electric  current  density  given 
by  expijno).  In  Fig.  4.  the  normalized  magnitude  of  Rn 
is  plotted  versus  n  for  a  1-wavelength  diameter  cylinder. 
As  we  might  expect,  residuals  corresponding  to  the  zeroth-, 
first-,  and  second-order  harmonics  are  zero,  indicating  perfect 
absorption,  due  to  the  fact  that  metrons  given  by  (2)  are 
equivalent  to  these  three  harmonics.  However,  it  is  evident 
from  the  figure  that  higher  order  current  distributions  result  in 
much  larger  residuals.  For  the  sake  of  comparison,  Fig.  4  also 
shows  the  normalized  magnitude  of  Fourier  spectrum  of  the 
actual  current  distribution.  No  significant  current  components 
are  present  in  the  region  where  the  residual  attains  large 
values.  We  have,  therefore,  verified  that  for  small  cylinders, 
or  whenever  a  current  distribution  to  be  computed  can  be 
expressed  as  a  linear  combination  of  metrons,  the  MEI  method 
yields  excellent  results. 

But  what  happens  if  the  aforementioned  conditions  are  not 
met.  as  might  be  the  case  with  larger  cylinders?  Fig.  5  is  a  plot 
of  the  magnitude  of  the  surface  current  density  due  to  a  TE- 
polarized  plane  wave  incident  on  a  cylinder  10  wavelengths 
in  diameter.  The  MEI  and  series  solutions  are  compared,  and 
we  observe  that  the  error  is  much  more  noticeable  than  in 
Figure  3  for  the  1-wavelength  diameter  cylinder.  However, 


Fig.  3.  Plot  of  Iht  magnitude  of  the  surface  current  density  on  a  PEC  circular 
cylinder  as  a  function  of  6  (D  =  lA,  Ap  =  A/20.  AS  =  b.~°.  2  layers.  TE 
polarization)  The  MEI  solution  is  compared  to  the  series  solution. 


Fig  4.  Plot  of  the  normalized  value  for  [Hnl  and  plot  of  the  normalized 
magnitude  of  the  Fourier  specUum  of  the  surface  current  density  for  the  case 
shown  in  Fig  3. 

the  MEI  solution  still  follows  the  general  shape  of  the  series 
solution. 

A  glance  at  the  spectrum  of  the  current  distribution  in  Fig.  6 
indicates  that  significant  harmonic  components  are  present  up 
to  approximately  the  35th  order.  The  metrons,  which  are  only 
associated  with  the  first  three  harmonics,  cannot  accurately 
represent  the  actual  current  on  the  cylinder,  yet  the  results  in 
Fig.  5  are  satisfactory.  Furthermore,  a  plot  of  the  normalized 
magnitude  of  the  residual  for  different  harmonics  in  Fig  6 
reveals  that  the  residual  is  practically  zero  for  all  the  harmonics 
up  to  approximately  the  20th  order  even  though  one  would 
expect  it  to  be  zero  only  for  the  first  three  harmonics  that  are 
used  as  metrons.  It  is  also  evident  that  the  residual  become-' 
significant  in  the  region  where  there  are  stiil  significant  specuai 
components  of  the  current.  This  explains  the  deterioration  in 
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Fig  5.  Plot  of  the  magnitude  of  the  surface  current  density  on  a  PEC  circular  *^'8  Flo'  of  magnitude  of  the  surface  current  density  on  a  PEC  circular 
cylinder  as  a  function  of  6  (D  =  ]0A,  Ap  =  A/20,  A0’=  0.5T°,  2  layers,  cylinder  as  a  function  of  6  (D  -  lOOX.Ap  =  A/20.  AW  =  0.0  .  2  layers. 
TE  polarization)  The  MEI  solution  is  compared  to  the  series  solution.  ^  polarization).  The  MEl  solution  is  compared  to  the  series  solution 


Fig.  6  Plot  of  the  normalized  value  for  |/?n  |  and  plot  of  the  normalized  ^*8  8  Flot  of  the  normalized  value  for  |2?,i  |  and  plot  of  the  normalized 
magnitude  of  the  Fourier  spectrum  of  the  surface  current  density  for  the  case  magnitude  of  Fourier  spectrum  of  surface  current  density  for  the  case  shou  n 
shown  in  Fig.  5.  Fig  7. 


accuracy  when  compared  to  the  smaller  problem  shown  in 
Fig.  3. 

The  overlap  between  regions  of  nonzero  residual  and 
nonzero  current  spectrum  becomes  even  more  pronounced 
as  the  size  of  the  problem  is  increased.  Fig.  7  is  a  plot  of  the 
magnitude  of  the  surface  current  for  a  cylinder  that  is  100 
wavelengths  in  diameter  due  to  a  TE-polarized  plane  wave.  In 
Fig.  8,  we  show  the  normalized  magnitudes  of  corresponding 
harmonic  residuals  and  current  spectrum. 

Let  us  summarize  the  most  important  observations  made 
thus  far.  First,  whenever  a  linear  combination  of  metrons 
accurately  approximates  the  actual  current  distribution,  the 
MEl  method  gives  excellent  results  since  the  residual  is 
zero  for  all  spectral  components  of  current  that  are  of  in¬ 
terest.  Second,  even  if  the  above-mentioned  conditions  are 


not  satisfied,  the  results  obtained  with  the  MEl  method  are 
very  good.  Although  the  residuals  for  only  the  firsi  three 
harmonics  arc  explicitly  forced  to  zero,  the  residual  is  rela¬ 
tively  small  for  most  of  the  remaining  harmonic  components 
of  interest.  In  other  words,  currents  that  are  much  more 
rapidly  varying  than  the  metrons  produce  fields  that  are  readil> 
absorbed  by  the  boundary  condition.  Third,  as  the  size  of 
the  problem  is  increased,  the  accuracy  is  impaired  ResiduaK 
corresponding  to  the  most  rapidly  varying  current  componenis 
increase. 

It  is  evident  by  now  that  the  accuracy  of  the  MEl  solu¬ 
tion  can  be  judged  by  observing  variations  in  the  harmonic 
residuals.  On  the  other  hand,  the  behavior  of  these  residuals 
is  essentially  determined  by  the  choice  of  metrons.  However, 
before  we  proceed  into  further  analysis  of  the.se  relationships. 
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Fij;  >).  Semili'gjnthmic  ploi  of  \P„  |  for  a  PEC  circular  cylinder  (D  =  lOA. 
A/i  =  .\/2U.  =  0.j7°,  2  layers.  TM  polarization). 


a  numerical  study  of  certain  theoretical  concepts  introduced  in 
Section  11  of  this  paper  is  considered. 

B  Suit  Field 

Throughout  the  remainder  of  this  paper  we  shall  be  referring 
to  results  of  the  MEI  method  as  applied  to  a  10-wavelength 
diameter  circular  cylinder  illuminated  with  a  TM-polarized 
plane  wave.  However,  this  is  merely  a  workhorse  for  our 
numerical  study  since  we  have  found  that  the  conclusions  hold 
for  circular  and  square  cylinders  of  various  sizes. 

Fig.  9  is  a  plot  of  the  magnitudes  of  the  harmonic  residuals 
when  the  low  order  sinusoidal  metrons  are  used  to  produce 
the  MEI  coefficients.  A  semilogarithmic  scale  reveals  certain 
details  not  \isible  in  the  conrespyonding  linear  scale  plot  shown 
by  the  solid  curse  in  Fig.  12.  Only  the  first  three  harmonics, 
that  is.  only  metrons.  yield  residuals  that  are  exactly  zero 
(finite  numbers  due  to  numerical  roundoff).  Three  regions 
in  the  variation  of  harmonic  residuals  are  now  evident.  The 
first  region,  where  the  residuals  are  exactly  zero,  coincides 
with  the  metrons.  Within  the  second  region,  the  residuals  are 
not  exactly  zero,  but  they  are  insignificant  for  all  practical 
purposes.  The  reason  for  the  existence  of  this  region  will  be 
explained  in  this  and  the  following  subsection.  The  fact  that 
the  MEI  produces  an  excellent  boundary  condition  for  the 
harmonics  in  this  region  is  the  reason  for  its  success  when 
applied  to  large  geometries.  Finally,  after  a  certain  threshold 
spatial  frequency,  the  residual  abruptly  increases,  indicating 
the  failure  of  the  boundary  equation  to  absorb  fields  radiated 
by  these  rapidly  varying  currents.  The  abrupt  increase  is  not 
discernible  on  the  semi-log  plot,  but  it  is  very  evident  on  the 
linear  plot  in  Fig.  12.  It  is  the  overlap  of  this  third  region  with 
the  region  of  nonzero  spectrum  of  the  current  distribution  that 
ultimately  determines  the  amount  of  error  in  the  MEI  solution. 

Based  on  the  above  discussion  of  the  behavior  of  the 
harmonic  residuals,  we  can  provide  some  additional  insight 
into  the  theoretical  concept  of  the  null  field,  defined  by  (7). 
The  null  field  is  simply  a  linear  combination  of  individual 


Fig.  10.  Plol  of  (he  magnitudes  of  the  Green’s  functions  due  to  line  source 
located  at  the  six  boundary  nodes  and  plot  of  the  magnitude  of  the  null  held 
for  the  case  shown  in  Fig.  9. 

fields  radiated  by  unit  line  sources  placed  at  all  the  nodes  that 
are  coupled  through  the  MEI  equation.  This  situation  is  shown 
in  Fig.  10.  The  magnitudes  of  all  six  line  source  fields  are 
plotted  along  with  the  magnitude  of  the  null  field  that  has  been 
obtained  by  linearly  combining  these  six  fields  through  the  use 
of  the  MEI  coefficients.  Fields  are  sampled  along  the  cylinder 
surface.  According  to  (14),  the  null  field  should  preferably  be 
as  small  as  possible  in  order  for  the  MEI  method  to  be  accurate. 
Indeed,  as  we  can  see  in  Fig.  10.  the  null  field  is  much  smaller 
than  any  of  the  six  components.  It  is  the  purpose  of  the  MEI 
coefficients  to  weight  each  one  of  the  line  source  fields  in  such 
a  way  that  they  cancel  each  other  as  throughly  as  possible. 

The  null  field  bears  a  very  close  relationship  to  the  harmoniv 
residuals.  In  fact,  since  R„  is  the  residual  due  to  nth  harmonic 
current,  it  is  easy  to  show  from  (12)  that  the  Rn  is  the  same 
as  the  nth  term  in  the  Fourier  scries  of  the  null  field.  It  i' 
therefore  possible  to  obtain  the  harmonic  residuals  simply  by 
finding  the  Fourier  spectrum  of  the  null  field.  Follow  ing  th.^ 
direction,  we  show  in  Fig.  1 1  the  magnitude  of  the  Founer 
spectra  for  the  fields  radiated  by  line  sources  placed  at  the 
six  nodes  on  the  mesh  boundary  where  the  MEI  equations 
are  to  be  applied.  One  thing  is  immediately  brought  to  oar 
attention.  The  spectra  arc  almost  identical  up  to  approximately 
the  20th  order,  where  the  harmonic  residuals  were  found  :o 
be  significant.  Indeed,  we  argue  that  it  is  exactly  this  prope'iy 
of  the  Green's  functions  emerging  from  nearby  points,  i.e.. 
that  they  differ  significantly  only  in  higher  spatial  frequencies, 
which  is  responsible  for  the  high  accuracy  of  the  MEI  equation 
even  when  applied  to  currents  that  are  much  more  rapidly 
varying  than  the  metrons  and  are  not  expressible  as  the  linear 
combination  of  the  latter.  We  shall  make  this  point  clearer  in 
the  next  section. 

C.  Parameter  Study 

In  the  MEI  method,  the  coefficients  were  chosen  such  thai 
the  first  three  harmonics  of  the  null  field  are  zero.  Howe\er. 
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I  ig  II.  Ploi  of  the  magnitudes  of  the  Fourier  spectra  of  six  Green's  functions 
mat  were  shown  in  Fig.  10. 


since  the  spectra  of  the  six  Green’s  functions  are  virtually 
identical  up  to  the  20th  harmonic,  we  should  be  able  to  use 
any  arbitrary  set  of  harmonics  between  0  and  20  and  obtain 
essentially  the  same  coefficients.  In  Fig.  12  we  see  that  this  is 
actually  the  case.  The  null  field  spectra  are  plotted  for  different 
choices  of  metrons.  The  traditional  choice  of  metrons,  i.e., 
zeroth,  first,  and  second  harmonic  as  metrons,  yields  almost 
an  identical  spectrum  to  that  obtained  from  the  use*  of  the 
zeroth,  fifth,  and  tenth  harmonics  as  metrons.  On  the  other 
hand,  the  use  of  metrons  outside  the  interval  where  Green’s 
function  spectra  are  identical,  as  for  example  using  the  zeroth, 
24th.  and  35th  harmonic  as  metrons,  gives  radically  different 
results  that  are  actually  better. 

For  the  sake  of  comparison,  Fig.  13  shows  the  magnitude 
of  error  in  the  MEl  solution  for  the  above-mentioned  meuon 
choices.  It  is  evident  that  the  solution  is  insensitive  to  the 
choice  of  metrons  as  long  as  they  are  chosen  from  the  spatial 
frequency  interval  where  the  Green’s  function  spectra  do  not 
differ  much.  It  is  also  evident  that  a  smaller  residual  leads  to 
a  better  solution. 

One  obvious  way  to  reduce  residual  and  thereby  improve 
the  solution  is  to  move  the  boundary  of  the  mesh  farther 
away  from  the  cylinder.  One  can  perform  a  simple  analysis  to 
show  that  the  differences  between  the  Green’s  functions  due 
to  the  six  line  sources  decrease  as  the  boundary  of  the  mesh 
is  moved  farther  from  the  cylinder.  Evidence  of  this  is  shown 
in  Fig.  14,  where  the  magnitudes  of  harmonic  residuals  are 
plotted  for  the  mesh  consisting  of  two,  four,  and  six  layers. 
The  correspionding  magnitude  of  the  error  in  the  current  on 
the  cylinder  is  shown  in  Fig.  15.  Thus,  moving  the  boundary 
farther  from  the  scatterer  reduces  the  residual  and  improves 
the  solution. 

There  are  other  ways  to  improve  the  solution,  such  as 
selecting  a  suitable  set  of  metrons.  We  have  already  seen  in 
Fig.  13  that  by  choosing  the  zeroth,  24th.  and  35th  harmonics 
to  .be  the  metrons,  we  obtain  better  results  than  for  the  other 
two  choices  that  were  investigated.  One  way  to  systemize  this 
choice  is  to  use  more  metrons  than  unknown  coefficients — in 


Fig.  12.  Plot  of  |/?n|  for  the  case  shown  in  Fig.  9.  but  for  three  different 
sets  of  metrons  ^roth,  first,  and  second  harmonics;  then  the  zeroth,  fifth, 
and  tenth  harmonics;  and  finall>  the  zeroth.  24th.  and  3Sth  harmonics  were 
used  as  metrons. 


Fig  13.  Plot  of  the  magnitude  of  the  error  in  the  MEI  solution  for  the  surface 
current  density  as  a  function  of  0  for  the  three  cases  show  n  in  Fig  12 


this  case,  more  than  five — and  then  to  solve  for  the  coefficients 
with  the  least  square  approach.  Fig.  16  shows  what  happens 
with  the  magnitude  of  harmonic  residuals  as  the  number  of 
meu-ons  is  increased  from  5  to  50  and  100.  In  Fig.  17.  we 
demonstrate  the  effect  of  this  upon  the  magnitude  of  error  in 
the  MEI  solution.  It  decreases  with  an  increasing  number  of 
metrons,  though  the  rate  of  decrease  in  the  solution  is  not  as 
fast  as  the  rate  of  decrease  in  the  residual. 

It  is  natural  now  to  ask  what  happens  if  we  keep  increasing 
the  number  of  metrons  to  infinity.  Would  the  MEI  coefficients 
be  stable  and.  if  they  are.  what  would  they  converge  to?  Using 
an  infinite  number  of  metrons  intuitively  corresponds  to  the 
requirement  that  the  residual  given  by  (13)  be  as  small  as 
possible  for  all  possible  currents.  That  brings  us  back  to  the 
argument  expressed  in  (14).  whose  solution  was  found  to  be 
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Fij  14.  Plot  of  |f?„|  for  the  case  shown  in  Fig  9.  but  for  three  different 
separations  between  the  boundar>  of  the  mesh  and  the  cylinder  surface.  Results 
for  two.  four,  and  six  layers  of  nodes  are  shown. 


Fig  15.  Plot  of  the  magnitude  of  the  error  in  the  MEl  solution  for  the  surface 
current  density  as  a  function  of  d,  for  the  three  cases  shown  in  Fig.  14. 


(15).  In  other  words,  one  could  use  five  metrons  equal  to  the 
conjugates  of  the  Green’s  functions  to  obtain  the  same  solution 
as  if  one  were  to  use  an  infinite  number  of  sinusoidal  metrons. 
This  is  verified  numerically  in  Table  1. 

Fig.  18  shows  the  magnitudes  of  the  harmonic  residuals 
corresponding  to  3{X)  sinusoidal.  5(X)  sinusoidal,  and  5g* 
metrons.  The  convergence  towards  the  g’  solution  is  evident. 

The  sinusoidal  metrons  usually  perform  either  as  well  as 
or  better  than  the  5*  metrons  for  the  cases  shown  here. 
We  expect  this  behavior  because  the  problem  of  plane  wave 
scattering  from  circular  cylinders  results  in  surface  currents 
whose  spectral  content  is  located  at  the  lower  harmonics.  Since 
the  g'  metrons  optimize  the  coefficients  for  all  the  harmonics, 
the  residual  for  the  lower  order  harmonics  are  not  as  well 
minimized  as  for  the  case  where  the  low-order  sinusoidal 
metrons  are  used.  It  should  be  noted,  however,  that  for  sources 
and  cylinder  shapes  that  induce  surface  currents  with  a  broader 
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Fig  16.  Plot  of  |/?„  I  for  the  case  shown  in  Fig.  9.  but  using  5.  50.  and  100 
lowest  order  harmonics  as  metrons. 


Fig.  17.  Plot  of  the  magnitude  of  the  error  in  the  MEI  solution  for  the  vurface 
cuireni  density  as  a  function  of  d  for  the  three  cas^-  shown  in  Figure  16. 


TABLE  I 


Type  of  metrons 

Number  of  metrons 

MEI  coefficient  03 

sinusoidal 

5 

-0.4821 -FjO.8103 

sinusoidal 

500 

-0.3027-)-j0.0520 

9' 

5 

-0.2993-Fj0.0464 

range  of  harmonics,  we  expect  the  g’  metrons  to  perfonn  better 
than  the  sinusoidal  metrons.  The  choice  of  metrons  is  really 
dependent  on  the  geometry  and  excitation  in  the  problem. 

IV.  Conclusion 

In  this  paper,  a  theoretical  foundation  has  been  pres  ided 
for  the  MEI.  It  was  shown  that  the  method  is  a  local  approx- 


1105 


JEVnC  AND  LEE  THE  MEASURED  EQUATION  OF  INVARIANCE 


otiler  ofhiiiDOnic 


Fig  18.  Plot  of  |/?„  |for  the  cast  shown  in  Fig.  9,  but  using  300  hannonics, 
500  hannonics.  and  the  conjugates  of  Green's  functions  as  meirons. 

* 

imate  boundary  condition.  However,  because  this  boundary 
condition  incorporates  geometrical  information,  it  can  be  ex¬ 
pected  to  produce  better  results  than  other  local  boundary 
conditions.  A  numerical  study  has  been  performed  showing  the 
influence  that  various  parameters  have  on  the  MEI  solution. 
The  parameters  include  the  electrical  size  of  the  scatterer, 
the  distance  between  the  mesh  boundary  and  the  scatterer, 
and  the  choice  of  metrons.  The  theoretical  analysis  has  been 
verified  by  the  numerical  results.  We  have  shown  why  the  MEI 
method  successfully  handles  current  distributions  that  are  not 
expressible  as  linear  combination  of  metrons.  We  have  also 
identified  the  major  sources  of  error.  It  was  shown  that  the 
performance  of  the  boundary  condition  can  be  predicted  before 
the  actual  problem  is  solved.  We  have  also  proposed  a  new  set 
of  metrons  that  gives  the  same  results  as  those  obtained  using 
an  infinite  number  of  sinusoidal  metrons. 

The  results  presented  here  and  in  previous  papers  are  only 
the  initial  work  on  this  new  method.  Based  on  our  expieriences 
with  this  method,  we  believe  that  major  improvements  on  the 
method  can  be  made  in  the  future. 
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